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EXECUTIVE  SUMMARY 


A  wireless  sensor  network  (WSN)  is  an  autonomous,  self-organizing  network 
without  any  pre-established  infrastructure  or  centralized  administration.  WSNs  have  been 
used  for  a  wide  range  of  applications  where,  often,  the  main  goal  is  to  monitor  a  specified 
phenomenon.  One  important  monitoring  task  which  has  recently  caught  the  attention  of 
WSN  researchers  is  that  of  locating  a  signal  source  by  extracting  the  infonnation 
contained  in  the  received  signal.  Two  primary  approaches  to  source  localization  have 
evolved:  direction  of  arrival  (DOA)  based  techniques  and  time  difference  of  arrival 
(TDOA)  based  techniques 

This  thesis  proposes  a  least  squares  error  estimator  for  DOA  localization  which  is 
unbiased  when  the  noise  is  Gaussian-distributed  with  zero  mean.  This  estimator  solves  an 
over  determined  Vandennonde  system  of  equations  which  is  known  to  be 
computationally  efficient  and  accurate. 

Based  on  this  least  squares  error  estimator,  this  thesis  proposes  a  passive  source 
localization  scheme  which  exploits  the  non-line-of-sight  (NLOS)  signals  from  non- 
cooperative  sources.  The  proposed  solution  is  a  hybrid  DOA/TDOA  source  localization 
scheme  and  is  comprised  of  three  parts:  a  DOA  estimator,  an  association  algorithm  for 
the  identified  signal  bearings,  and  the  source  localization  scheme  itself.  The  recently 
proposed  Space  Division  Multiple  Access  (SDMA)-based  receiver  is  used  for  DOA 
estimation.  TDOA  information  is  used  to  discriminate  between  the  line-of-sight  (LOS) 
and  the  NLOS  signals  and  to  associate  the  incoming  multipath  signal  with  the 
corresponding  source  and  reflector  pair.  In  special  cases,  the  proposed  scheme  is  capable 
of  solving  the  association  problem  spatially  without  the  need  for  TDOA  information.  A 
technique  is  also  provided  to  estimate  the  position  and  the  orientation  of  the  reflectors 
when  site-specific  database  information  is  not  available.  Both  centralized  and  distributed 
variants  of  the  proposed  scheme  are  presented  with  the  latter  being  of  particular  interest 
in  WSNs. 


The  proposed  localization  scheme  allows  a  wireless  sensor  network  to  (1)  perform 
single  array  localization,  (2)  perform  the  localization  in  a  distributed  fashion,  (3)  obtain 
source  location  estimates  with  NLOS  signals  only  and  (4)  improve  the  location  estimates 
compared  with  those  obtained  using  the  LOS  infonnation  only. 

A  simulation  model  was  developed  and  implemented  in  MATLAB  to  evaluate  the 
performance  of  the  proposed  localization  scheme.  The  simulation  results  demonstrate  that 
the  proposed  localization  scheme  provides  high  accuracy  estimates  and  outperfonns  the 
LOS-only  based  localization  schemes.  This  is  primarily  because  more  bearings  are 
available  and  the  conditioning  of  the  least  squares  problem  is  better  for  the  proposed 
scheme.  Furthermore,  the  simulation  results  also  show  that  the  proposed  scheme  is  able 
to  provide  accurate  NLOS-only  source  location  estimates  when  the  LOS  paths  are  not 
available. 
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I.  INTRODUCTION 


A.  INTRODUCTION  TO  PASSIVE  SOURCE  LOCALIZATION  USING 

WIRELESS  SENSOR  NETWORKS 

A  wireless  sensor  network  (WSN)  is  an  autonomous,  self-organizing  network 
without  any  pre-established  infrastructure  or  centralized  administration  [1].  WSNs  have 
been  used  for  a  wide  range  of  applications  where,  often,  the  main  goal  is  to  monitor  a 
specified  phenomenon  [2].  Wireless  sensor  networks  offer  numerous  advantages  when 
compared  to  traditional  wired  or  wireless  networks  [3].  Of  particular  interest,  WSNs 
provide  greater  redundancy  since  the  malfunction  of  a  number  of  sensors  has  less  impact 
on  the  overall  system  perfonnance.  Additionally,  WSNs  can  be  deployed  quickly  at  low 
cost  and  are  well-suited  to  use  in  mobile  platforms.  Not  surprisingly,  they  have  found 
wide-spread  interest  in  emergency  and  military  applications. 

One  important  monitoring  task  which  has  recently  caught  the  attention  of  WSN 
researchers  is  that  of  locating  a  signal  source  by  extracting  the  information  contained  in 
the  received  signal.  That  source  could  be  an  enemy’s  radio  location,  as  in  military 
applications,  or  the  direction  of  arrival  (DOA)  and  the  location  of  a  sensed  phenomenon 
in  an  emergency  situation  like  the  seismic  waves  which  follow  an  earthquake  [3].  A 
WSN  performs  the  source  localization  by  coordinating  the  effort  of  individual  sensors 
which  act  as  antenna  elements.  These  sensors  are  clustered  together  to  fonn  antenna 
arrays  which  fuse  the  data  collected  by  the  sensors,  to  carry  out  the  localization  task  [4]. 

B.  RELATED  WORK  IN  PASSIVE  SOURCE  LOCALIZATION  WITH 

WIRELESS  SENSORS  NETWORKS 

Passive  source  localization  using  wireless  sensor  arrays  is  a  problem  addressed 

extensively  in  existing  literature.  Two  primary  approaches  to  source  localization  have 

evolved  [5]:  direction  of  arrival  (DOA)  based  techniques  and  time  difference  of  arrival 

(TDOA)  based  techniques.  DOA-based  localization  systems  utilize  antenna  arrays  which 

examine  the  spatial  characteristics  of  the  incident  signals  to  obtain  bearing  estimates  [11], 

[12],  [13],  [14],  [15]  and  [16].  The  bearing  estimates  are  used  for  the  position  location 

determination  by  triangulation  [5]  [17],  [18]  and  [19].  TDOA  localization  systems 

1 


estimate  the  source  location  using  the  intersection  of  hyperboloids,  which  are  the  set  of 
range  difference  measurements  between  three  or  more  receiving  sensors.  These  are 
determined  by  measurement  of  the  TDOA  of  a  signal  between  the  sensors  [2],  [6],  [7], 
[8],  [9],  and  [10]. 

The  need  for  accurate  localization  requires  knowledge  of  the  spatial 
characteristics  of  the  wireless  channel  since  those  characteristics  significantly  affect  the 
performance  of  the  arrays.  Thus,  a  significant  challenge  is  the  development  of  realistic 
channel  models  which  can  accurately  predict  the  behavior  of  the  wireless  medium  [30]. 
In  a  real-world  deployment,  complex  propagation  phenomena  lead  to  uncertainty  in 
deciding  whether  a  direction  of  arrival  (DOA)  corresponds  to  the  line-of-sight  (LOS) 
signal  or  its  reflection  [2],  That  uncertainty  can  lead  to  significant  errors  when  estimating 
the  position  of  the  desired  source.  The  problem  is  more  severe  when  a  LOS  signal  does 
not  exist  as  is  often  the  case  in  urban  environments.  Additionally,  it  has  been  reported 
that  a  number  of  strong  reflections  can  be  expected  even  in  rural  areas  [31],  These  non- 
line-of-sight  (NLOS)  signals  potentially  provide  additional  information  that  can  be 
exploited  in  the  source  location. 

There  have  been  several  proposals  in  literature  that  consider  the  presence  of 
NLOS  signals.  The  first  category  is  comprised  of  schemes  which  attempt  to  mitigate  the 
effects  of  the  NLOS  signals.  In  [22],  the  measurements  are  weighted  to  emphasize  the 
LOS  signals,  while  [23]  identifies  the  arrays  that  do  not  receive  LOS  signals  and  excludes 
them  from  the  localization  process.  In  both  approaches,  the  authors  seek  to  minimize  the 
impact  of  the  NLOS  signals  rather  than  taking  advantage  of  them.  Recently,  a  second 
category  of  solutions  is  beginning  to  emerge  that  attempts  to  exploit  these  NLOS  signals. 
The  authors  of  [24]  propose  a  hybrid  DOA/TDOA  scheme  which  exploits  the  NLOS 
when  the  desired  source  is  cooperative. 

In  contrast,  this  thesis  proposes  a  passive  source  localization  scheme  which 
exploits  the  NLOS  signals  from  non-cooperative  sources.  It  is  a  DOA/TDOA-based 
scheme  which  uses  bearing  estimation  for  localization  through  triangulation.  The  TDOA 
information  is  used  to  discriminate  between  the  LOS  and  the  NLOS  signals  and  to 


2 


associate  the  incoming  multipath  signal  with  the  corresponding  source  and  reflector  pair. 
Furthermore,  using  the  NLOS  information,  the  proposed  scheme  is  capable  of  performing 
single  array  localization. 

C.  THESIS  OBJECTIVE 

The  objective  of  this  research  is  to  develop  a  source  localization  scheme  that  is 
capable  of  non-cooperative  source  localization  within  the  constraints  of  a  WSN.  The 
existing  work  on  cooperative  source  localization  can  then  be  viewed  as  a  special  case  of 
this  more  general  solution. 

There  are  two  significant  contributions  in  this  work.  The  first  is  a  proposed  least 
squares  estimator  for  DOA-based  localization.  The  second  is  a  passive  source  localization 
scheme  which  exploits  the  NLOS  signals  from  non-cooperative  sources.  The  latter  is  a 
DOA/TDOA-based  scheme  which  uses  bearing  estimation  for  localization  through 
triangulation.  The  TDOA  information  is  used  to  discriminate  between  the  LOS  and  the 
NLOS  signals  and  to  associate  the  incoming  multipath  signal  with  the  corresponding 
source  and  reflector  pair.  Furthermore,  using  the  NLOS  information,  the  proposed 
scheme  is  capable  of  performing  single  array  localization  or  NLOS  only  based 
localization.  Both  centralized  and  distributed  variants  of  the  proposed  scheme  are 
presented  with  the  latter  being  of  particular  interest  in  WSNs. 

D.  THESIS  OUTLINE 

Chapter  II  provides  the  background  to  support  the  proposed  work.  It  introduces 
the  fundamental  concepts  of  antenna  arrays  including  the  response  of  an  array  with 
randomly  distributed  elements.  An  overview  of  the  wireless  environment  is  also  included 
to  validate  the  adopted  propagation  and  received  signal  models.  The  chapter  then  outlines 
several  current  approaches  to  DOA  estimation.  The  recently  proposed  Space  Division 
Multiple  Access  (SDMA)-based  receiver  is  presented  and  compared  to  the  conventional 
MUSIC  algorithm.  Chapter  II  closes  with  a  brief  discussion  of  TDOA  estimation. 

The  significant  contributions  of  this  work  are  presented  in  Chapter  III.  Examining 

LOS-only  localization  first,  a  least-squares  estimator  for  DOA  localization  is  proposed 

and  analysis  is  provided  to  investigate  the  biasing,  the  impact  of  errors,  and  the 
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conditioning  of  the  proposed  estimator.  A  comparison  between  the  least-squares  and  the 
total-least-squares  solutions  is  also  presented.  To  take  advantage  of  a  long  duration 
signal,  a  sequential  least-squares  approach  is  also  included.  Turning  our  attention  to 
NLOS  as  well  as  LOS  signals,  the  proposed  localization  scheme  which  exploits  the 
NLOS  signals  is  then  described.  The  incident  signal-source-reflector  association 
algorithm  is  outlined  and  a  technique  is  provided  to  estimate  the  position  and  the 
orientation  of  the  reflector.  Chapter  III  concludes  with  the  source  localization  procedure 
of  the  proposed  scheme. 

In  Chapter  IV,  simulation  results  are  provided  to  evaluate  the  performance  of  the 
proposed  scheme  and  compare  it  to  existing  LOS-only  solutions. 

Chapter  V  summarizes  the  significant  contributions  of  this  thesis  and  provides 
some  ideas  for  extending  this  work  in  the  future. 

Finally,  the  Appendix  includes  the  MATLAB  code  used  in  the  simulation. 


4 


II.  BACKGROUND  ON  SOURCE  LOCALIZATION 


This  chapter  provides  the  background  to  support  the  proposed  solution  for  passive 
source  localization  of  non-cooperative  sources  using  clusters  of  random  or  aperiodic 
sensor  arrays  which  form  a  wireless  sensor  network  (WSN).  In  this  thesis,  we  consider  a 
hierarchical  WSN  that  consists  of  two  levels.  The  top  level  is  the  network  of  arrays  in 
which  each  array  is  viewed  as  a  single  node.  These  nodes  perform  the  DOA  estimation 
task  and  the  bearings  obtained  are  transmitted  to  a  fusion  center,  where  the  localization 
algorithm  is  executed.  In  the  proposed  distributed  variant  of  our  solution,  each  node  also 
performs  single-array  localization.  In  this  case,  the  source  location  estimations  are  then 
simply  averaged  by  the  data  fusion  center.  In  the  second  level  of  the  WSN,  each  array 
contains  a  number  of  sensor  elements.  DOA  estimation  is  conducted  at  this  sensor  layer, 
while  the  source  localization  is  conducted  at  the  array  level.  An  example  WSN  is  shown 
in  Figure  1. 


Fusion  Center 


Source 


Figure  1 .  Example  WSN  comprised  of  3  arrays,  a  single  source  of  interest  and  a 

fusion  center. 
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This  chapter  begins  by  discussing  the  array  response  to  an  incident  signal.  An 
overview  of  the  propagation  environment  follows  and  a  propagation  prediction  model  is 
adopted  to  provide  a  realistic  scenario  for  performance  comparison  of  the  localization 
schemes.  A  link  budget  analysis  is  then  conducted  and  the  signal  received  by  the  array  is 
computed.  Finally,  this  chapter  closes  with  a  discussion  of  both  DOA  estimation  and 
TDOA  estimation.  Two  methods  of  DOA  estimation  are  presented,  the  Multiple  Signal 
Classification  (MUSIC)  DOA  estimation  method  [26]  and  the  new  Space  Division 
Multiple  Access  (SDMA)  receiver  [25].  A  comparison  of  the  angle  resolution  and  the 
accuracy  of  both  algorithms  is  included  and  the  SDMA  receiver  is  chosen  as  the  DOA 
estimation  method  of  the  localization  process  in  this  thesis.  The  TDOA  estimates  will  be 
used  to  associate  the  multipath  components  of  the  received  signal  with  the  corresponding 
source  and  the  related  reflectors. 


A.  ARRAY  RESPONSE 

The  array  steering  vector  contains  the  responses  of  all  sensors  to  a  source  with  a 
single  frequency  of  unit  power  [26].  The  array  response  varies  as  a  function  of  direction 
and  a  steering  vector  is  associated  with  each  direction  of  interest.  The  uniqueness  of  the 
association  is  defined  by  the  array  configuration  [27].  In  this  section,  we  begin  with  a 
discussion  of  linear  arrays  followed  by  the  more  general  case  of  two-dimensional 
aperiodic  and  random  arrays. 


1.  Uniform  Linear  Array  (ULA) 

The  configuration  of  a  uniform  linear  array  (ULA)  is  shown  in  Figure  2.  The 
source  transmits  a  narrow  band  signal  .s'  ( t )  of  frequency  /  and  is  assumed  to  be  in  the  far 

field.  The  array  consists  of  N  sensors  with  unifonn  inter-element  spacing  d  .  With 
respect  to  the  reference  node  (sensor  1),  sensor  2  experiences  a  time  delay  of  [28] 

d  cos  6 


At  =  ■ 


(1) 


where  c  is  the  signal  propagation  speed.  The  time  delay  At  corresponds  to  a  phase  shift 
of  the  signal  equal  to 


6 


(2) 


=  2n 


d  cos  6 

I 


This  phase  shift  is  the  same  for  every  pair  of  sensors  because  the  inter-element  spacing  is 
constant  (uniform).  Assuming  identical  sensor  elements,  the  steering  vector  of  this  array 
is  given  by 


a(<9)=  1 


eiKv 


e~j(  N-l)A*r 


(3) 


which  can  in  turn  be  written  as 
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(«) 


i 


-  i—d  cos  6 

e  1 


-  j(N-\)—dcosd 

e  * 


(4) 


s(0 


Figure  2.  Uniform  linear  array  of  N  sensors  with  inter-element  spacing  d.  The 

source  is  located  in  the  far  field. 

2.  Two-Dimensional  Aperiodic  and  Random  Arrays 

The  geometry  of  a  two-dimensional  aperiodic  array  is  shown  in  Figure  3.  The 
sensors  of  the  array  are  placed  in  th ex-y  plane  according  to  some  algorithm.  Without  a 
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loss  of  generality,  the  position  of  the  reference  sensor  is  assumed  to  be  at  the  origin  of  the 
coordinate  system.  The  phase  difference  between  sensor  i  and  the  reference  sensor  is 
given  as 


2  n 

T 


(x;.  cos  6  +  y,  sin  6*) 


(5) 


and  thus  the  corresponding  steering  vector  is  found  to  be 


a 


(») 


i 


-j — (jc2  cos  0+y2  sin#) 

e  a 


2.7T 

cos  0+yt  sin  0 ) 


(6) 


When  the  positions  of  the  sensors  are  chosen  by  some  random  process,  the  aperiodic 
array  is  known  as  random  array  [29].  The  steering  vector  of  the  random  array  is  identical 


to  that  of  the  aperiodic  array  except  that  the  vector 


1  e 


2k 

cos0+y2  sin#) 


is  a  random  vector  as  in 


2  7i 

COS @+yi  s*n 


(7) 


Figure  3.  Two-dimensional  array  geometry. 


This  thesis  considers  aperiodic,  random  arrays.  Aperiodic  sensor  arrays  have 
several  advantages  when  compared  to  conventional  periodic  sensor  arrays.  Due  to  the 
non-periodic  nature  of  the  sensor  spacing,  they  do  not  suffer  from  grating  lobes  in  their 
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spectrum  and  are  not  limited  to  a  maximum  sensor  separation  of  0.5  X  (0.5  m  at  300 
MHz)  [29].  This  reduces  the  cost  of  the  construction  of  such  an  array,  since  a  smaller 
number  of  sensors  are  required.  Larger  element  separation  also  provides  more  resilience 
against  mutual  coupling  [29]  which  occurs  between  the  sensors  when  one  is  in  the 
vicinity  of  the  other  [26].  This  coupling  effect  degrades  array  performance  and  is  largely 
ignored  in  most  array  signal  processing  literature.  Finally,  random  arrays  provide 
flexibility  in  deployment  and  can  accommodate  arbitrary  topologies,  which  are  common 
in  mobile  platforms  and  WSNs. 

B.  SOURCE  TO  ARRAY 

In  this  section,  we  discuss  and  adopt  a  radio  propagation  model  of  the  wireless 
communications  channel  that  will  be  used  in  the  subsequent  performance.  Free  space 
path  loss  causes  signal  strength  attenuation  of  a  LOS  electromagnetic  wave,  while  with 
NLOS,  multipath  components  can  also  be  attenuated  due  to  reflection,  diffraction,  and 
scattering  [31].  Electromagnetic  signals  experience  attenuation  while  they  travel  in  space. 
This  is  the  result  of  spherical  energy  spread  in  space.  Reflection  occurs  when  a  signal 
strikes  a  surface  and  is  then  reflected  towards  the  receiver.  Diffraction  is  the  phenomenon 
that  occurs  when  the  electromagnetic  signal  strikes  the  edge  of  the  corner  of  a  large 
structure  compared  with  the  wavelength.  Scattering  occurs  when  a  signal  strikes  an  object 
that  is  much  smaller  compared  with  the  wavelength  [32],  These  effects  lead  to  complex 
multipath  propagation  scenarios,  especially  in  areas  with  a  high  density  of  potential 
reflectors  and  scatterers  (e.g.,  in  urban  areas). 

1,  Propagation  Environment 

In  this  section,  we  look  at  signal  attenuation  in  both  line-of-sight  and  reflected 
signals.  We  present  loss  expressions  for  both. 

a.  Free  Space  Loss 

In  wireless  communications,  as  the  signal  propagates  through  the  medium, 
it  disperses  with  distance  [33].  This  type  of  attenuation,  known  as  free  space  loss  ( Ls ) , 

can  be  expressed  as  the  ratio  of  the  signal  power  between  the  transmitter  and  the  receiver 
in  dB  as 
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where  Pt  is  the  transmitted  signal  power,  P  is  the  received  signal  power,  and  d  is  the 

distance  between  transmitter  and  receiver.  The  free  space  loss  is  proportional  to  the 
square  of  the  distance  between  the  transmitter  and  the  receiver.  Thus,  as  this  distance  is 
increased,  the  free  space  attenuation  becomes  very  large. 


b.  Reflection 

When  a  signal  propagating  in  one  medium  encounters  the  boundary  of 
another  medium,  it  is  partially  reflected  back  to  the  first  medium  and  partially  refracted  to 
the  new  medium  [21].  The  propagation  characteristics  of  the  resulting  waves  are 
governed  by  the  boundary  conditions  of  the  interface. 

A  schematic  representation  of  the  reflection  on  a  smooth  surface  is  shown 
in  Figure  4.  The  relationships  for  the  angles  shown  in  the  figure  are  given  by  Snell’s  law 
as 

9,  =  9„  and 

r-  (9) 

sin  0R  =  fsr  sin  0T 


where  sr  is  the  dielectric  constant  of  the  material.  The  attenuation  of  the  reflected  signal 
is  given  by  the  square  of  the  reflection  coefficient  as  in  PR  =  Pt  |r|" .  The  reflection 

coefficient,  T  ,  is  a  function  of  the  reflected  and  refracted  angles  and  has  a  range  between 
0  and  1.  This  coefficient  depends  on  the  type  of  polarization.  Thus,  for  transverse- 
electric  (TE)  polarized  plane,  it  can  be  shown  to  be  [21] 

cos  9„  -  cos  9r 

rE  =  *  V - (10) 

cos  0R  +  cos  0T 

while  for  transverse-  magnetic  (TM)  polarized  plane  it  is  [21] 

xfe~  cos  0„  -  cos  0T 

r  H=^= - - - (ii) 

f£r  cos  0R  +  cos  6t 
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For  er  =  5 ,  Figure  5.  plots  the  absolute  value  of  the  reflection  coefficient 
for  both  types  of  polarization  as  a  function  of  the  incident  angle.  For  TM  polarization, 
there  exists  a  single  incident  angle  at  which  TH  =  0  and  no  reflection  occurs.  This  angle 
is  called  the  Brewster  angle  and  is  expressed  as  [3 1] 


sin 


K)  = 


4Er~X 


(12) 


On  the  other  hand,  when  the  incident  angle  is  90°,  |r£|  =  |Tff  |  =  1  and  the 

signal  is  totally  reflected  for  both  types  of  polarization.  The  path  loss  because  of  the 
refection  is 


4=-20iog(|r|) 


(13) 


Figure  4.  Geometrical  representation  of  the  reflection.  The  incident  signal  with 
angle  0X  is  partially  reflected  with  angle  0R  and  partially  refracted  with 

angle  6r 
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Figure  5.  |F|  as  a  function  of  the  incident  angle  for  sr  =  5  . 

2.  Radio  Propagation  Prediction 

Shifting  our  attention  to  the  multipath  effects,  the  classic  radio  propagation 
models  provide  information  about  the  signal  power  and  Doppler  shifts  of  the  received 
signal  [30].  However,  in  source  localization  schemes,  the  time  delay  spread  and  the 
angle-of-arrival  spread  are  also  of  major  importance.  A  fundamental  challenge  in 
localization  applications  is  therefore  the  development  of  realistic  channel  models  that  can 
accurately  predict  the  characteristics  of  the  multipath  propagation.  These  propagation 
models  are  highly  dependent  upon  the  propagation  environment  and  there  does  not  exist  a 
single  model  that  covers  every  environment.  Rather,  several  empirical  models  have  been 
proposed  based  on  measurements  from  different  areas  within  the  specific  environment  of 
interest.  In  general,  the  characterization  of  these  environments  is  based  upon  the 
population  density  and  the  building  architecture  (e.g.,  urban,  suburban,  and  rural).  The 
most  challenging  propagation  environment,  particularly  for  source  localization  problems, 
is  the  urban  environment  where  the  received  signal  may  be  dominated  by  strong 
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reflections  and  the  LOS  path  is  not  always  present.  A  typical  propagation  scenario  can  be 
seen  in  Figure  6  and  this  thesis  primarily  focuses  on  this  source  localization  problem  in 
urban  environments. 


Figure  6.  Typical  propagation  scenario.  A  LOS  path  and  several  NLOS  paths  are 

shown. 

a.  Site-Specific  Propagation  Prediction 

Most  propagation  prediction  models  for  urban  areas  examine  the  statistical 
parameters  of  the  propagation  environment,  such  as  average  row  spacing  and  building 
height  distributions  [21].  That  methodology  is  suitable  when  the  aim  of  the  model  is  the 
prediction  of  the  average  power  of  the  signal  received.  However,  applicability  is  limited 
in  localization  applications  where  the  localization  scheme  is  actively  attempting  to 
exploit  the  NLOS  components  of  the  signal. 

As  an  alternative,  the  site-specific  prediction  model  takes  advantage  of  the 
actual  mapping  of  the  area  under  consideration.  The  mapping  can  be  exploited  by  using  a 
database  which  contains  the  footprints  of  the  buildings  and  any  other  large  objects  [21]. 
Those  databases  must  be,  in  general,  three-dimensional.  In  the  special  case  when  the 
assumption  of  low  antennas  in  tall  building  environments  holds,  the  databases  can  be 

13 


two-dimensional.  This  is  because,  in  those  cases,  the  primary  propagation  paths  lie 
around  the  sides  of  the  buildings  and  no  significant  propagation  paths  exist  over  the 
rooftops  [21]. 

Site-specific  propagation  models  use  ray  optical  methods  and  treat  each 
reflection  as  an  individual  ray  in  the  space.  The  maximum  number  of  reflections  along  a 
ray  path  could  be  anywhere  from  six  to  eight  and  the  number  of  walls  in  the  database  on 
the  order  of  thousands  [21]  —  facts  which  make  the  prediction  a  complex  procedure. 
Additionally  though  (theoretically)  the  total  received  power  is  given  by  the  summation  of 
the  powers  of  the  individual  rays;  this  strategy  does  not  give  a  very  accurate  estimate 
because  it  ignores  the  scattered  signals.  Finally,  the  accuracy  of  the  databases  is  on  the 
order  of  0.5  m  in  the  best  case  [21],  which  introduces  further  uncertainty  into  the 
prediction. 

All  of  this  notwithstanding,  a  site-specific  approach  offers  a  simplifying 
feature  when  applied  to  source  localization.  When  a  ray  path  includes  many  reflections, 
the  power  delivered  by  the  ray  is  reduced  dramatically  due  to  the  attenuation  of  each 
reflection  and  the  increased  path  length.  In  many  cases,  only  the  first  reflections  are 
strong  enough  to  be  exploited.  Thus,  the  power  of  each  individual  ray  is  more  important 
in  the  DOA  estimation  problem  than  the  total  received  power  and  ignoring  the 
contribution  of  the  weak  components  can  be  an  advantage.  The  localization  procedure  is 
significantly  simplified  by  ignoring  the  multi-reflected  rays,  which  could  not  be  exploited 
anyway. 

With  respect  to  the  size  of  the  reflector,  there  is  a  limitation  on  the  size  of 
the  region  over  which  reflections  occur  and  ray  methods  can  be  used.  This  can  be 
understood  by  looking  at  Figure  7  in  which  a  ray  from  the  source  is  reflected  by  the 
reflector  as  shown.  The  length  of  the  reflector  surface  is  wB ,  while  2 wF  is  the  width  of 
the  Fresnel  zone,  defined  by  the  relative  geometry  between  the  receiver  and  the  reflector. 
If  wB  cos  6  >  2  wF ,  then  the  rules  for  reflection  can  be  used  in  order  to  predict  the 

propagation  characteristics  [21],  However,  since  max  {2 wF}  =  VXs  where  5*  is  the  total 
length  of  the  reflected  ray  in  meters,  the  maximum  ray  length  is  then  bounded  by 
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(14) 


S  <  —(wB  cos  6)" . 

A 

Thus,  if  9  -  45° ,  wB  =  20  m ,  and  A  =  0.5  m  ,  (14)  limits  the  ray  length  to  400  m  .  If  the 
above  inequality  does  not  hold,  then  the  object  acts  as  a  scatterer  [34]. 


Figure  7.  Fresnel  zone  and  reflector  surface  characteristics  [From  21]. 

b.  Scattering  Prediction 

The  propagation  rays,  as  discussed  in  the  previous  section,  were 
considered  discrete  lines  in  space  which  possess  the  total  power  of  the  corresponding 
multipath  component.  However,  in  reality,  the  power  of  the  signal  is  distributed  around 
these  rays  which  deteriorate  the  accuracy  of  the  DO  A  estimation  methods.  The 
distribution  of  the  power  around  the  ray  is  the  result  of  the  signal  scattering  caused  by 
small  objects  found  around  the  source.  The  same  phenomenon  is  observed  around  the 
reflections  which  can  be  considered  secondary  sources.  Thus,  the  scatterers  are  grouped 
into  clusters,  around  both  the  source  and  the  reflectors. 

Several  models  have  been  proposed  in  literature  to  predict  this  scattering. 
The  distribution  of  the  angles  of  arrival  caused  by  the  scattering  can  be  modeled  as 
discrete  Gaussian  [30]  or  discrete  Laplacian  [35]  as  in 
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Gaussian 


(15) 


e 2 

PG{0)  =  cle2*2 

/  \  -'/2- 

PL  [0)  =  c2e  CT  Laplacian 

where  6  <=\-7t,+n}  and  the  mean  value  of  both  functions  is  the  angle  which  corresponds 

to  the  associated  ray.  The  parameter  a  controls  the  spread  of  the  functions.  Typically,  a  is 
small  and  the  values  of  both  functions  are  concentrated  around  the  mean.  A  plot  of  the 
function  for  a  =  3  is  shown  in  Figure  8.  Measurements  in  [36]  indicate  that  the 
distribution  of  the  incident  angle  fits  the  discrete  Laplacian  function  better,  since  in  both 
rural  and  urban  environments,  they  tend  to  demonstrate  a  sharp  peak  while  also 
maintaining  long  tails. 


Figure  8.  Gaussian  versus  Laplacian  scattering  functions. 

3.  The  Adopted  Propagation  Model 

In  summary,  the  propagation  model  used  in  this  thesis  combines  elements  from 
the  majority  of  the  models  discussed  above.  In  a  site-specific  propagation  model,  the 
power  attenuation  of  each  ray  is  calculated  as  the  sum  of  the  path  loss  and  the  reflection 
attenuation.  The  site-specific  mapping  provides  the  number  of  reflectors  encountered  by 
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the  ray  en  route  to  its  destination  and  the  path  loss  is  computed  using  the  free  space  path 
loss  model  for  each  ray  individually.  The  total  power  of  the  received  signal  is  the 
summation  of  the  powers  of  the  individual  rays. 

C.  RECEIVED  SIGNAL 

The  received  signal  model  assumes  multiple,  uncorrelated  sources  transmitting 
signals  to  an  N-sensor  array.  In  this  section,  a  link  budget  analysis  is  performed  and  the 
received  signal-to-noise  ratio  (SNR)  of  the  incoming  signals  is  computed.  Additionally,  a 
matrix  expression  is  defined  for  the  received  signal  as  a  function  of  the  response  of  the 
receiving  array. 

1.  Link  Budget  Analysis 

The  power  level  of  the  incident  signals  is  computed  by  subtracting  the  losses  due 
to  attenuation,  as  described  previously,  from  the  transmitted  power  and  by  adding  the 
gain  of  the  array.  The  gain  of  the  array  is  given  in  dB  as  [29] 

G„=101og(W)  +  G£I  (16) 

where  GEL  is  the  gain  of  each  array  element  in  dB.  Expression  (16)  holds  for  both  ULA 

and  aperiodic  arrays,  although  in  the  latter  case,  it  must  be  considered  an  approximation 
[37].  The  link  budget  equation  in  dB  for  each  signal  path  is  then  expressed  as 

SNR aj,  =  SNRsc  -Ls-Lr+  Gar  +Gt  (17) 

where  SNR4R  is  the  signal-to-noise  ratio  available  at  the  receiving  array  for  a  specific 
signal  path,  SNRsc  is  the  signal-to-noise  ratio  of  the  signal  transmitted  from  the 
corresponding  source  Ls  is  the  free  space  loss,  LR  is  the  path  loss  because  of  the 
reflection,  GT  is  the  transmitter  gain  (assumed  to  be  0  dB  in  this  work)  and  G4R  is  the 
array  gain.  The  SNR  for  signal  power  P  is  defined  in  dB  as 
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SNR  =  10  log  - 

\kTB 


(18) 


where  k  is  the  Boltzmann  constant  equal  to  1.38x10  23  nr  kg  s~2  K1 ,  T  is  the  system 
noise  temperature  in  degrees  Kelvin  which  includes  both  the  antenna  and  receiver  noise 
and  B  is  the  effective  noise  bandwidth  of  the  receiver. 

2.  Received  Signal  Model 

Although  the  signal  of  each  source  is  considered  narrowband  in  the  previous 
discussion,  the  results  can  be  extended  to  wideband  signals,  given  the  assumption  that  the 
frequency  response  of  the  array  is  flat  over  the  signal’s  bandwidth  and  the  propagation 
time  across  the  array  is  small  when  compared  to  the  inverse  of  the  bandwidth  [26].  If  the 
number  of  sources  is  K  and  each  of  them  transmits  a  signals*  (t) ,  the  received  signal  at 
time  t  in  the  array  can  be  expressed  as 

*M=Z“KhM+"W  09) 

k= 1 

where 

J(t)  =  [xj(t)  x2  (?)  ...  e  CiV ,  (20) 

n  (t)  is  the  noise  vector  and  a(0k)  is  the  steering  vector  for  signal  sk  (t)  with  DOA  6k . 
In  matrix  notation,  (20)  can  be  written  as 

x(t)  =  {t)  +  n  (0  (21) 

where  A (O  j  e  C V/A  is  the  array  response  matrix  which  contains  the  responses  of  each 
sensor  for  each  incident  signal  and  is  given  as 

A[0)  =  [a{0x)  a{02)  ...  a(0k)\  (22) 

and  6  is  a  vector  which  contains  the  DOA  for  each  incoming  signal.  Let  L  be  the 
number  of  observations  with  L>  K  .  The  received  signal  can  then  be  written  as 
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A  =  .4(0)S  +  A, 


(23) 


where  N  e  CN/L  is  the  noise  matrix  and  S  e  CKxL  is  the  transmitted  signal  matrix  . 

D.  DOA  ESTIMATION 

This  section  discusses  and  compares  two  methods  of  DOA  estimation,  the 
Multiple  Signal  Classification  (MUSIC)  estimation  and  the  Space  Division  Multiple 
Access  (SDMA)  receiver. 

1.  The  MUSIC  Algorithm 


Most  DOA  algorithms  are  based  on  the  computation  of  the  signal  correlation 
matrix  [32],  The  received  signal  correlation  matrix  Rxx  and  the  desired  signal  correlation 

matrix  Rss  are  defined  as 


K=E{x(t)x-(t)} 

(24) 

«.  =  £{s(r)s"(«)} 

(25) 

where  H  denotes  the  Hermitian  transpose  of  the  matrix  and  E  ( a  j 

is  the  expectation  of 

the  argument  a  [26].  If  the  statistics  of  the  signal  and  the  noise  are  not  known  but  the 

corresponding  processes  are  ergodic,  then  the  correlation  matrices  can  be  approximated 

by  averaging  a  finite  number  of  data  observations  as  [26] 

4=7^x(f/K(0 

L  /= 1 

(26) 

^  l=\ 

(27) 

Both  Rss  and  Rxx  are  NxN  matrices.  If  additive  white  Gaussian  noise  is  assumed,  the  two 

matrices  are  related  as  [26] 

A=4®)kX(®)+<W 

(28) 

where  cr2  n  is  the  noise  variance  and  /  is  an  NxN  identity  matrix.  If  the  incident  signals 

are  uncorrelated,  Rss  is  diagonal.  When  the  signals  are  coherent,  Rss  is  singular  [38].  In 
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most  cases,  the  signals  are  partially  correlated  and  Rss  is  positive  definite.  This  property  is 
very  important,  since,  as  will  be  shown  later,  the  DOA  estimation  algorithms  are  based 
on  the  inversion  of  Rss.  Rxx  has  N  eigenvalues  (/fy  A,,  ...  , AN )  and  N  associated 


eigenvectors  E  =  [e1  e2  ...  eN\  which  can  be  obtained  by  the  eigenvalue 

decomposition.  By  ordering  the  eigenvalues  from  larger  to  smaller,  the  eigenvector 
matrix  can  be  divided  into  two  sub-matrices  [32] 


£  = 


E  E 


(29) 


These  sub-matrices  are  also  called  subspaces.  Es  has  K  columns  and  corresponds  to  the 

signal  subspace,  while  En  has  N -K  columns  and  corresponds  to  the  noise  subspace.  An 
alternate  way  to  find  the  eigenvectors  of  the  autocorrelation  matrix  is  directly  to  use  the 
received  signal  matrix  X  and  eigendecompose  it  by  using  the  singular  value 
decomposition  (SVD),  which  gives  more  stable  algorithms  [39]. 


The  Multiple  Signal  Classification  (MUSIC)  [26]  is  the  most  popular  among  the 
DOA  algorithms  based  in  the  subspace  decomposition  of  the  correlation  matrix.  The 
desired  DO  As  are  estimated  by  identifying  the  peaks  of  the  MUSIC  spatial 
pseudospectrum,  which  is  given  as  [40] 


,  aH  (d)a(d) 

“  =  aH{0)EnEHna{ey 


(30) 


The  array  steering  vectors  are  orthogonal  to  the  noise  subspace  and,  therefore,  the  peaks 
in  the  pseudospectrum  represent  the  DOAs  for  the  desired  signal. 

MUSIC  can  be  applied  to  any  arbitrary  but  known  array  topology  [26]  and 
accordingly,  requires  accurate  array  calibration  [17].  Furthermore,  MUSIC  assumes  the 
number  of  sources  is  known  in  order  to  assign  the  corresponding  eigenvectors  to  the 
signal  subspace  and  several  algorithms  have  been  proposed  to  do  this  [41].  Additionally, 
conventional  MUSIC  breaks  down  under  near-coherent  signal  conditions  like  those 
which  exist  when  multipath  propagation  conditions  are  present  [32].  Again,  several 
methods  have  been  proposed  to  address  this,  but  they  are  either  applicable  only  to  special 
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cases  such  as  spatial  smoothing  which  works  for  ULA  [17]  or  they  are  computationally 
intensive  such  as  multidimensional  MUSIC  [17]. 

2.  The  SDMA  Receiver 

The  SDMA  receiver  is  a  new  method  for  DOA  estimation  proposed  in  [25].  A 
depiction  of  this  receiver  is  shown  in  Figure  9.  The  SDMA  receiver  does  not  rely  on  the 
subspace  decomposition  of  the  correlation  matrix,  rather  it  cross-correlates  the  received 
signal  with  a  pre-computed  set  of  array  responses  for  every  direction  of  interest. 


Figure  9.  SDMA  receiver  [From  32], 


From  Figure  9,  consider  an  N  sensor  array.  The  output  of  each  array  element 
xn  (?)  is  phase-modulated  by  a  set  of  uncorrelated  spreading  sequences  wn  (7) .  These 

sequences  can  be  any  type  of  pseudorandom  [32]  or  orthogonal  [37]  sequence.  The 
produced  array  outputs  are  then  orthogonal  or  nearly  orthogonal.  In  matrix  notation,  the 
output  of  the  array  is  found  to  be 

Y  =  WhX  (31) 

where  X  is  given  from  (22)  and  W  e  C  Vx/'  is  the  matrix  of  the  spreading  sequences 
which  is  defined  as 

W  =  [wx(t)  v v2(t)  ...  wN(t)]H .  (32) 
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The  signals  stored  in  the  virtual  array  are  also  modulated  by  the  same  set  of  spreading 
sequences.  In  matrix  notation,  the  output  of  the  virtual  array  is  given  as 

V  =  WHAM(&),  (33) 

where  AM  (0)  e  CNxL  is  the  matrix  of  the  array  responses  of  all  sensors  for  all  DOAs.  The 
correlator  cross-correlates  the  array  signal  with  that  of  the  output  of  the  virtual  array  as  in 

R  =  VhY  (34) 

where  R  is  a  KxL  matrix.  The  spatial  spectrum  of  the  SDMA  receiver  is  then 

Psdma  =  |^a-|  (3^) 

where 

Rk=  2>(U)  1>(2,/)  ...  2>(^0  •  (36) 

_  /=1  /= 1  /= 1 

The  peaks  of  PSDMA  correspond  to  the  DOAs  of  the  incident  signals. 

The  array  used  by  the  SDMA  receiver  is  preferably  a  two-dimensional  random 
array  such  that  the  sensor  geometry  and  element  phasing  is  unique  for  each  DOA  [32]. 
The  spreading  technique  used  further  ensures  that  all  received  directions  of  interest  are 
uniquely  defined.  The  receiver  does  not  compute  the  correlation  matrix,  but  just 
correlates  the  received  signal  with  a  pre-computed  one.  Thus,  in  contrast  to  MUSIC,  it 
does  not  rely  on  the  correlation  matrix  properties.  Also,  the  SDMA  receiver  does  not 
require  knowledge  of  the  number  of  incident  signals.  Finally,  the  SDMA  receiver  does 
not  rely  on  any  complex  adaptive  or  slow  iterative  methods.  It  “looks”  in  pre-detennined, 
finite  set  of  directions  of  interest  [32]  which,  by  estimating  multiple  angles 
simultaneously,  translates  to  an  instantaneous  search  through  a  bank  of  a  finite  number  of 
expected  observations. 

3.  Comparison  between  SDMA  Receivers  and  MUSIC 

A  simulation  was  conducted  using  MATLAB  and  a  comparison  was  made  of  the 
resolution  and  the  accuracy  between  the  SDMA  receiver  and  the  MUSIC  algorithm. 
Figure  10  shows  the  correlation  magnitude  of  SDMA  and  MUSIC  for  seven  incoming 
signals  with  angular  separations  of  20  degrees.  A  random  array  of  31  sensors,  occupying 
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a  25  m2  rectangular  area,  was  used.  All  the  signals  had  a  carrier  frequency  of 
300  MHz  (l  =  lm)  with  a  SNR  of  15  dB .  The  results  were  averaged  across  100  Monte 
Carlo  simulations.  It  is  clearly  evident  that  even  though  the  number  of  signals  is  large, 
SDMA  provides  steep  lobes  in  the  directions  of  the  incident  signals,  while  MUSIC  gives 
main  lobes  which  are  difficult  to  identify.  Furthermore,  SDMA  slightly  outperforms 
MUSIC  in  terms  of  accuracy  as  well.  Table  1  provides  the  measured  mean  value  //  and 

the  variance  a1  in  degrees  of  the  average  error  of  two  incident  angles.  Again  the  results 
were  averaged  across  100  Monte  Carlo  simulations.  As  expected,  both  schemes  provide 
unbiased  estimates  since  the  mean  values  are  near  zero,  but  SDMA  has  a  slightly  smaller 
variance. 


Figure  10.  Performance  comparison  of  SDMA  and  MUSIC  for  random  array  of  3 1 
sensors  that  covers  an  area  of  25  m2 .  SNR  of  incident  signals  is  15  dB . 


A 

C72 

MUSIC 

0.0917 

0.6325 

SDMA 

-0.015 

0.5908 

Table  1 .  Error  performance  of  SDMA  and  MUSIC. 
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Using  the  MATLAB  histfit.m  command,  the  distribution  of  the  error  of  both 
schemes  was  compared  with  the  theoretical  zero-mean  Gaussian  distribution.  The  results 
are  shown  in  Figure  11  for  the  SDMA  receiver  and  in  Figure  12  for  the  MUSIC 
algorithm.  The  results  from  both  schemes  can  be  seen  to  be  closely  approximated  by  the 
error  distribution  which  is  indicated  by  the  red  curve  in  both  Figures  1 1  and  12. 


Figure  1 1 .  SDMA  error  distribution. 
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300 


Figure  12.  MUSIC  error  distribution. 


Finally,  Figure  13  shows  the  performance  of  the  SDMA  receiver  and  MUSIC  in 
the  case  of  sparse  arrays.  The  same  parameters  from  Figure  10  were  used,  except  that  the 
array  area  was  increased  to  2500  m2 .  A  significant  improvement  in  the  resolution  of  both 
methods  is  observed.  As  discussed  earlier,  this  is  an  important  advantage  when  using 
random  arrays.  The  larger  spacing  between  the  array  elements  corresponds  to  a  larger 
array  aperture,  which  provides  higher  angular  resolution  [42].  In  the  sparse-array  case, 
SDMA  is  also  seen  to  outperform  MUSIC. 


25 


1 


CD 

TD 

=3 

'E 

CD 

CD 


C 

o 

to 

0) 

o 

O 


—  SDMA  Recei\*r 

-  MUSIC 

- 

~ 

w 

‘■'■'A.'vVf- 

N- 

ft 

VM 

-80  -60  -40  -20  0  20  40  60  80 


DOA,  0 


Figure  13.  Performance  comparison  of  SDMA  and  MUSIC  for  a  random  array  of  3 1 
sensors  that  covers  an  area  of  2500  nr  .  SNR  of  incident  signals  is  15  dB  . 


E.  TDOA  ESTIMATION 

In  TDOA-based  localization  schemes,  the  first  step  estimates  the  time  difference 
of  arrival  (TDOA)  of  the  incident  signal  between  the  sensor  nodes  in  the  network.  TDOA 
information  can  be  obtained  by  two  general  methods:  the  first  involves  the  subtraction  of 
the  time  of  arrival  between  sensors  to  obtain  the  relative  difference  while  the  second 
method  uses  cross-correlation  techniques  to  estimate  the  desired  TDOA  [17].  The  first 
method  requires  knowledge  of  the  transmission  time  which,  in  the  case  of  non- 
cooperative  sources,  is  not  available.  Accordingly,  only  the  cross-correlation  method  will 
be  discussed  here. 

Assume  that  the  signal,  s(t ) ,  is  transmitted  by  an  unknown  source.  Each  sensor  in 
the  cluster  will  receive  amplitude-scaled  and  time-delayed  versions  of  the  transmitted 
signal  corrupted  by  the  noise  introduced  in  the  channel.  Hence,  the  received  signals  at 
two  different  sensors,  x{  (t)  and  x2  (t)  ,  are  given  by 
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(37) 


xl  (V)  =  Als{t-r^)  +  nl  (V)  and 
x2  (7)  =  A2s  (7  -  r2)  +  n2  (7) 

where  A1  and  A2  are  the  amplitude  scaling  factors,  nt(t)  and  n2(t)  are  the  additive 
noise,  and  r,  and  r7  are  the  signal  offset  times  at  each  sensor.  At  and  A2  reside  in  the 
interval  [0,l] .  Assumes  that  the  noise  is  zero-mean  Gaussian  and  the  signal  and  the  noise 
are  uncorrelated.  Equation  (37)  can  be  rewritten  in  the  form 

Xi(f)  =  s(f)  +  «i(f) 

x2  (?)  =  As  (7  -  r)  +  n2  (7) 

where,  without  loss  of  generality,  it  is  assumed  that  the  first  sensor  is  the  one  with  the 
smaller  time  of  arrival,  the  amplitudes  are  normalized  by  Al  and  r, ,  again  without  loss  of 

generality,  is  set  to  0.  The  cross-correlation  between  the  signals xx(t)  and  x2{t)  can  be 
approximated  by  the  estimate 

Rxy{r)  =  ^-\xl(t)x2{t-T)dt  (39) 

1  0 

where  T  is  the  time  interval  during  which  the  observation  is  conducted.  The  lag,  r  ,  which 
maximizes  (38)  is  the  estimate  of  the  TDOA  value  [43]. 

F.  SUMMARY 

This  chapter  provided  the  background  to  support  the  proposed  solution  for  passive 
source  localization  of  non-cooperative  sources.  It  began  by  discussing  the  array  response 
to  an  incident  signal.  An  overview  of  the  propagation  environment  followed  and  a 
propagation  prediction  model  was  adopted  to  provide  a  realistic  scenario  for  performance 
comparison  of  the  localization  schemes.  It  was  shown  that  through  a  link  budget  analysis, 
the  signal  received  by  the  array  can  be  computed.  Finally,  the  chapter  concluded  with  a 
discussion  of  both  DOA  and  TDOA  estimation  techniques. 

In  the  next  chapter  we  propose  a  least  squares  estimator  for  DOA-based 
localization.  Based  on  this  estimator  we  develop  a  hybrid  DOA/TDOA  localization 
scheme  which  exploits  both  the  LOS  and  the  NLOS  signals. 
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III.  PASSIVE  SOURCE  LOCALIZATION  USING  RANDOM 

SENSOR  ARRAYS 

In  this  chapter,  we  address  the  problem  of  passive  source  localization  and  present 
the  proposed  non-cooperative  source  localization  scheme.  In  the  context  of  random 
arrays,  we  begin  with  an  analysis  of  LOS-only  DOA-based  localization  and  then  move  on 
to  localization  using  both  LOS  and  NLOS  signals. 

A.  LOS-ONLY  DOA-BASED  LOCALIZATION 

DOA-based  localization  schemes  use  DOA  estimates  of  the  incident  signals  to 
obtain  bearings  to  the  position  of  a  source  of  interest.  The  bearing  estimates  are  used  for 
the  position  location  determination  by  triangulation.  In  a  two-dimensional  case,  there  are 
two  unknowns  (the  coordinates  of  the  source  position).  In  practice,  though,  more  than 
two  bearings  are  required  due  to  finite  angular  resolution,  multipath  fading,  and  noise 
[17].  These  bearings  fonn  an  over-determined  system  of  equations  which  often  do  not 
intersect  at  a  single  point  and  can  be  solved  using  the  least-squares  approach.  Several 
proposals  exist  in  literature  based  on  the  Least  Squared  Error  (LSE)  minimization  [18], 
[45],  while  others  involve  the  Minimum  Mean  Square  Error  estimation  or  the  Maximum 
Likelihood  [46]. 

1.  The  Geometry  of  the  Problem 

The  proposed  least  squares  estimator  for  DOA  localization  is  based  on  the 
geometric  configuration  of  Figure  14.  Consider  K  arrays  where  each  array  obtains  a  DOA 
estimate  by  processing  the  incoming  signal  from  the  source.  9k  is  the  bearing  obtained  by 
the  array  k.  We  begin  by  formulating  the  least  squares  solution  of  this  problem. 
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Figure  14.  LOS  DOA-based  source  localization  using  three  arrays. 


From  Figure  14, 


tan  9k 


y-yk 

x-xk  ’ 


which  can  be  rewritten  as 

y-x tan 6k  =  yk  - xk  tan 6k . 
Doing  this  for  all  the  arrays  yields,  in  matrix  form, 


1 

1  1 
p-  p- 

&  ^ 
_ i 

y 

y1  -  tan  6yx l 
y2  -  tan  02x2 

1  -  tan  0k 

X 

_}’k  ~  tan  9kxk_ 

(40) 


(41) 


(42) 


or 


A(0)z=b(O). 


(43) 


30 


2.  Least  Squares  Solution 

The  matrix  A  ( 0 )  in  (42)  and  (43)  has  a  specific  structure  and  is  referred  as  a 

Vandermonde  matrix.  The  solution  of  least  squares  problems  involving  such  a  matrix  is 
known  to  be  computationally  efficient  and  accurate  [47]. 

a.  Least-Squares  Estimator 

An  estimate  of  the  solution  to  (42)  is  obtained  by  minimizing  the  L2  -norm 
of  the  residual  [48] 

r  =b(0)-A(6)z  .  (44) 

This  is  graphically  depicted  in  Figure  15  where  y  is  the  point  on  the  range  of  A  which 
is  closest  to  b  .  From  the  nonnal  equations,  the  solution  to  (42)  can  be  obtained  as 

z=A+(0)b{0 )  (45) 

where  A+ (0)  =  \^AH  (O')  A(0)^  A(0)H  is  the  pseudo-inverse  matrix. 


Figure  15.  Graphical  representation  of  the  least  squares  problem  (From  [48]). 
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b.  Statistical  Analysis  of  the  Least-Squares  Estimator 

To  show  that  this  estimator  is  unbiased,  let  us  assume,  as  discussed 


previously,  that  the  DOA  measurements  contain  a  zero-mean  Gaussian  error  80  such  that 


0esl=0+80. 

(46) 

Correspondingly, 

4,(e)  =  A(0)  +  SA(0) 

(47) 

and 

ba(e)=b(e)+sb(e). 

(48) 

Taking  the  derivative  of  both  A(0 )  and  b  (0)  from  (42), 

dA(0) 

dO 


1 

xi 

cos2  0X 

cos2  0X 

1 

a l.  (  n\ 

X2 

cos2  02 

at)  (u) 

and  - = 

cos2  02 

dO 

1 

xk 

cos2  0k 

cos2  6k 

~  dO,  8b  a 

db  and  8  A  a  dA  , 

5  A  = 


0 

r 

i 

80xxx 

cos2  0X 

cos2  0X 

0 

S02 

S02x2 

cos2  01 

and  8b  = 

cos2  02 

0 

sok 

80kXk 

cos2  0k 

cos2  0k 

Substituting  (47)  and  (48)  into  (43) 

\Tf~r  r-T~ 


( A  +  8A)t  (b  +  8b)  =  (A  +  8A)r  (A  +  8A)(z+8z ) 


(49) 


(50) 


(51) 


where  8z  is  the  error  in  the  source  location  estimation.  Expanding  (51)  and  ignoring  the 
second  order  terms, 
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ATb  +  AT 5b  +  5ATb  =  AT Az  +  AT5Az  +  5ATAz  +  AT  AS  z  . 
Rearranging  terms, 

Ar(5b- 5Az)  +  [AT  +5 AT)[b  -Az)  =  AT AS z  . 
With  b  -  Az  =  0,  reduces  to 

(ATA)lAT(5b-5Az)  =  5z. 

Using  the  pseudoinverse  we  finally  arrive  at 

5z  =A+(5b-5Az) 


which  has  the  expectation 


E{5z}  =  A+E{5b-5Az 


The  vector  inside  the  expectation  of  (53)  is 


5b  -  5Az  = 


( x  -  x, ) 

V  ,  x}  59x 
cos  9X 

M®, 

cos'  02 


(■ ) 

cos2  6, 


50, 


and  has  a  mean  value  of 


E 1 5b  -  5Az 


cos2  01 

(x-X2) 
COS2  On 


(■ x~xk ) 

cos2  0, 


E{S6 ;} 
E{592] 


E{56k 


(52) 


(53) 


(54) 


(55) 


However,  E{86k}  =  0  for  every  since  the  DOA  error  is  zero- 

mean  and,  therefore,  E^Sb  -5Azj  =  0.  From  (53),  it  follows  that  E  (<5Tj  =  0  and  the 
proposed  estimator  is  unbiased. 
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The  accuracy  of  the  proposed  scheme,  and  of  that  of  any  DOA-based 
localization  scheme,  is  affected  by  the  error  in  DOA  estimation.  A  metric  of  the  accuracy 
of  these  schemes  can  be  explained  in  tenns  of  the  covariance  matrix  [49] 

C,=(A«Ci]Ay  (56) 

where  Cs  e  CKxK  is  the  covariance  matrix  of  the  DOA  error  estimates  (assumed  to  be 
zero  mean  white  noise  sequence  in  this  case).  Cf)  is  a  diagonal  matrix  and  each  element 

along  the  main  diagonal  represents  the  error  variance  of  the  measured  DOA  for  the 
corresponding  array.  If  the  statistics  of  the  error  sequence  are  known  and  each  array  has  a 
different  noise  variance,  appropriate  weighting  of  the  least-squares  solution  leads  to  an 
optimum  estimator  [49].  The  weighting  matrix  is  chosen  to  be  the  C0 1  and  the  general 
form  of  the  weighted  least-squares  estimator  is 

xa=(A"C,'Ay  A"Cs'b.  (57) 

c.  Conditioning  of  the  Least-Squares  Solution 

The  conditioning  of  the  least-squares  problem  captures  the  perturbation 
behavior  of  the  least-squares  solution  [47].  A  least-squares  problem  is  characterized  as 
either  well-conditioned  or  ill-conditioned.  An  ill-conditioned  problem  is  one  in  which  a 
small  perturbation  in  the  observed  data  leads  to  a  large  error  in  the  estimated  solution  and 
the  condition  number  is  a  metric  used  to  quantify  the  conditioning.  A  well-conditioned 
problem  has  a  condition  number  equal  to  one.  As  the  condition  number  increases,  the 
problem  becomes  increasingly  ill-conditioned. 

In  the  least-squares  problem  of  (42),  the  condition  numbers  which 
describe  the  solution  z  ,  with  respect  to  the  perturbations  of  the  matrix  A  and  the  vector 
b  as  shown  in  Figure  15,  are  given  by  the  equations  of  Table  2  [47]  where 

k^A)  =  —,  l</r(yl)<oo  , 
cr2 

■  v  n 

9  =  cos  A=-  ,  0  <  9  <  —  , 
b  2 
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A  z 

r|  =  ,  0  <ij<k 

HI 

and  cr,  and  a2  are  the  maximum  and  the  minimum  singular  values  of  A ,  respectively. 

Here,  the  rank  of  the  matrix  A  is  2  and,  therefore,  two  singular  values  exist.  The 
conditioning  of  the  least-squares  problem  clearly  depends  on  the  geometry  of  the  specific 
scenario  under  consideration. 


Condition  number 

b 

77  cost? 

A 

,  .  k{A\  tan# 
k(A)  +  — - — - - 

77 

Table  2.  Condition  numbers  of  the  least-squares  problem  [From  47]. 

d.  Geometric  Dilution  of  Precision 

The  geometry  fonned  by  the  arrays  which  participate  in  the  localization 
procedure  and  the  source  of  interest  affects  the  accuracy  of  the  estimated  source  position 
[5].  As  the  distance  between  the  source  and  the  baseline  increases  relative  to  the  length  of 
the  baseline,  the  accuracy  of  the  least-squares  solution  decreases.  This  is  because  larger 
relative  distances  correspond  to  smaller  differences  in  bearings  at  the  arrays.  Thus, 
relatively  large  source  distances  increase  the  condition  number  of  matrix  A  and  tend  to 
make  it  singular.  This,  in  turn,  affects  the  stability  of  the  least-squares  problem.  This 
phenomenon  is  called  geometric  dilution  of  precision  and  is  illustrated  in  Figure  16.  As 
the  relative  distance  of  the  source  from  the  baseline  is  increased,  the  uncertainty  area, 
formed  by  the  intersections  of  the  estimated  bearings  is  also  increased.  Furthermore,  the 
angle  between  the  orthogonal  bisection  of  the  baseline  and  the  bearing  to  the  source  from 
the  midpoint  of  the  baseline  also  affects  the  accuracy  of  the  solution.  As  illustrated  in 
Figure  17,  as  angle  0BS  increases,  the  uncertainty  area  becomes  larger. 
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Figure  16. 


Figure  17. 


Array  1 


Array  2 


The  effect  of  distance  to  the  source  on  the  geometric  dilution  of  the 

precision. 


/ 


Array2 


Array  1 


Array2 


The  effect  of  the  bearing  to  the  source  on  geometric  dilution  of  the 

precision. 
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3. 


Total  Least  Squares 


In  the  proposed  least  squares  estimator,  both  A  and  b  can  be  affected  by  errors. 
For  a  least-squares  problem,  the  total  least  squares  (TLS)  method  is  known  to  compensate 
for  errors  in  matrix  A  [49].  Let  C  =  [  4  h  J  with  C  e  Ctx3  .  The  singular  value 
decomposition  of  matrix  C  is 

C  =  UZVH  (58) 


where  (7  e  Clx3,  E  e  M3x3  and  f  e  C3x3 .  Matrices  U  and  V  are  unitary  and  Z  is  a 
diagonal  matrix  of  the  form 


0 

°2C 

0 


0 

0 

^3  c 


(59) 


where  its  diagonal  elements  are  the  singular  values  of  C ,  satisfying  oxc  >  o2C  >  cr3C . 
The  TLS  solution  of  the  least-squares  problem  using  the  total  least-squares  estimator  can 
then  be  shown  to  be  [49] 

*TLS  =  (AH  A  ~  ^C1)  1  A"b  ■  (60) 

A  simulation  was  conducted  in  MATLAB  to  compare  the  performance  of  the 
proposed  model  using  both  the  LS  and  the  TLS  methods.  The  DOA  estimation  error  was 
assumed  to  be  zero-mean  Gaussian  with  a  standard  deviation  of  0.5  degrees.  The  results 
are  shown  in  Figure  18  as  a  function  of  the  number  of  sensor  arrays.  The  LS  solution  can 
be  seen  to  outperform  the  TLS  solution  for  small  numbers  of  arrays,  while  the  TLS 
solution  is  slightly  better  for  10  arrays  or  more. 
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Figure  18.  RMS  error  for  LS  vs.  TLS  for  1000  Monte  Carlo  simulations  with 

°doa  =°-5°- 

These  results  are  surprising  and  contradict  existing  literature.  As  mentioned  in 
[49],  the  TLS  is  a  deregularizing  procedure.  The  conditioning  of  the  TLS  problem  should 
always  be  worse  than  the  conditioning  of  the  respective  LS  problem.  The  ratio  of  the 
smallest  singular  value  of  A  to  the  smallest  singular  value  of  C  has  been  defined  as  a 
metric  to  quantify  the  instability  of  the  TLS  and  to  indicate  whether  the  LS  outperforms 

the  TLS  [49].  As  this  ratio,  in  this  case,  approaches  unity,  the  TLS  tends  to  be 


unstable  and  the  LS  solution  performs  better.  This  ratio  was  computed  for  the  simulation 
of  Figure  18  and  is  shown  in  Table  3.  The  results  can  be  seen  to  agree  with  the 
performance  of  the  TLS  algorithm. 


Arrays 

3 

4 

5 

6 

7 

8 

9 

10 

11 

a2A 

5.4513 

5.4934 

6.3222 

6.9740 

7.8332 

8.9594 

9.9635 

11.0359 

12.0155 

°'iC 

Table  3.  Stability  measure  for  the  TLS  algorithm. 
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4.  Sequential  Least  Squares 


As  time  progresses,  more  incoming  signal  samples  may  arrive  at  the  array  from 
the  same  source.  Hence,  more  data  will  be  available  and,  by  exploiting  it,  the  estimate  of 
the  source  location  can  be  improved.  In  this  case  the  localization  procedure  becomes  a 
dynamic  task  and  an  algorithm  which  continuously  updates  the  information  matrices  is 
needed  [48].  If  the  noise  is  uncorrelated  (which  implies  that  the  covariance  matrix  C()  is 


diagonal),  then  A  can  be  computed  sequentially  as 

'  A[k- 1] 


A[k\  = 


aH[k ] 


(61) 


where  A\k\  is  the  data  matrix  of  k measurements,  A [k  -l]  is  the  k  x2  matrix  of  the 


previous  k  - 1  measurements  and  a"  [/c]  is  the  kth  measurement.  Recall  that 


— H 

a 


[k]  =  [l  tan  6k  ]  for  the  proposed  scheme.  Thus,  the  sequential  estimator  is  [48] 


z[k]  =  z[k-l]  +  D[k](b[k]-aH[k]J[k-l]) 


where 


A  *]  =  ■ 


Cg  [k-l]a  [k] 


c j2k  +aH  [k] Cg  [k - 1] a  \k\ 
and  the  covariance  update  is  given  as 

Cg[k]  =  (l-D[k]a»[k])Cg[k-l]. 


(62) 


(63) 


(64) 


B.  LOS  AND  NLOS  DOA-BASED  LOCALIZATION 

The  proposed  hybrid  DOA/TDOA  LOS  and  NLOS  localization  scheme  is 
implemented  in  three  steps.  A  block  diagram  of  the  proposed  scheme  is  shown  in  Figure 
19.  The  first  step  implements  the  DOA  estimation  using  the  SDMA  receiver  as  described 
in  the  previous  chapter.  A  detailed  description  of  steps  2  and  3  (bearing  association  and 
source  localization,  respectively)  are  included  in  the  following  subsections. 
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Step  1 


Step  2 


Step  3 


Figure  19.  Block  diagram  of  the  proposed  LOS  and  NLOS  localization  scheme. 

1.  Association  of  Bearings 

Once  the  DOAs  have  been  determined  for  the  incident  signals,  they  are  then 
associated  with  each  other  according  to  source-reflector  pairs.  Specifically,  for  each 
source-reflector  pair,  a  LOS  and  a  NLOS  is  identified  based  on  knowledge  of  the 
reflector  position  and  geometry.  This  reflector  knowledge  can  be  obtained  either  from 
static  databases  which  contain  the  footprint  of  the  large  objects  within  the  environment 
under  consideration  [21]  or  dynamically  by  sending  a  beacon  from  a  known  location  and 
using  the  reflector  position  estimation  algorithm  which  is  described  in  the  following 
subsection.  When  neither  of  these  tools  is  available,  the  proposed  scheme  is  capable  of 
performing  reflector  position  estimation.  In  this  case,  the  problem  becomes  one  in  which 
the  scheme  simultaneously  performs  both  reflector  mapping  and  source  localization.  In 
the  remainder  of  this  section,  we  begin  with  a  discussion  of  how  this  reflector  mapping  is 
accomplished  and  then  present  the  association  scheme. 

a.  Reflector  Position  Estimation 

To  discuss  the  reflector  mapping  algorithm,  consider  Figure  20  where  a 
single  source  transmits  a  narrowband  signal  and  N  reflectors  generate  NLOS  signals  that 
can  be  viewed  as  secondary  sources.  A  minimum  of  three  arrays  must  be  available  to 
solve  the  reflector  estimation  problem.  Let  0.  n  be  the  bearing  associated  with  array  i  and 
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reflector  n  and  0ls  the  LOS  bearing  from  source  S  to  array  i.  The  bearing 


correspondences  can  be  illustrated  with  a  [N  +  4,3N  +  3) -bipartite  graph  as  shown  in 

Figure  21.  There  are  N  +  4  vertices  in  the  graph  and  3N  +  3  edges.  The  number  of  edges 
corresponds  to  the  total  number  of  bearings.  The  interrelationship  between  the  bearings  is 
determined  in  two  steps  as  shown  in  Figure  22.  In  the  first  step,  all  the  intersections 

between  the  bearings  of  array  1  and  array  2  are  found.  In  this  step,  a  total  of  (iV  +  1)” 

2x2  systems  of  equations  are  solved.  In  step  two,  the  residuals  between  those  points  and 
the  bearings  of  array  3  are  computed.  The  N  + 1  smaller  residuals  indicate  which 
intersections  of  array  1  and  array  2  bearings  correspond  to  which  bearing  of  array  3.  This 

second  step  involves  the  solution  of  (iV  +  1)'  equations.  Thus,  the  total  number  of 


required  operations  grows  as  N7 .  The  source  and  the  image  sources  related  to  the 
reflectors  are  estimated  by  solving  for  the  interrelated  bearings  using  the  least  squares 
estimator  of  the  proposed  scheme.  The  orientation  and  the  position  of  reflector  n 
according  to  Figure  20  is  then  given  as 


with 


z„  = 


y+y 


Rn 


2 

x  +  x'Rn 


0.  = 


x  =  x 


Rn 


arctan 


y  Rn~y, 

KX\n-Xsj 


-90° 


Rn 


(65) 


(66) 
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Figure  20. 


Figure  21. 


Unknown  reflector  position  and  orientation  estimation  with  three  arrays. 


Representation  with  a  (N  +  4,3N  +  3)  biparte  graph  of  the  bearing 
correspondences  of  the  reflector  source  pairs  for  N  reflectors. 
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Figure  22.  Two-step  procedure  for  finding  the  bearing  correspondences  of  the 

reflector  source  pairs 


b.  Bearing  Association  using  Expected  TDOA  Estimation 

Once  the  positions  of  the  reflectors  have  been  detennined,  each  array  can 
associate  its  bearings  using  a  set  of  time  markers.  From  Figure  23  and  using  the  sinusoids 
law, 

dlNC  _  _ ^ REF _  _  _ dLQS _  (67) 

sin (0*  -  0, )  sin (0R  +  0L- 20 o )  sin {20 R  - 20 o ) ' 

where  dAA,  is  the  distance  between  the  array  and  its  image  (with  respect  to  the  reflector), 

0O  is  the  orientation  of  the  reflector,  0,  is  the  LOS  bearing  and  0R  is  the  NLOS  bearing. 

However, 

Ad  —  dLOS  —  dINC  —  d  Rh:F  (68) 

and,  combining  (67)  and  (68), 
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(69) 


Ad 


sin  (20 R  -  20 0 )  -  sin  (dR  +  dL  -  20 0 )  -  sin {0R  -  0L )  ^ 
&m(0R+0L-20o) 


We  also  note  that 


d 


REF 


'AA' 


2s  in(0R-0o) 


and 


(70) 


where  c  is  the  propagation  speed  and  A  t  is  the  desired  expected  TDOA  between  the  LOS 
and  the  NLOS  signals.  Substituting  (69)  into  (70)  the  expression  for  At  is 

<P 


At  = 


sin  [0r+6l  -2do^m[0R  ~0o) 


dAA'^ 

2  c 


where 


(p  =  sin(2^ 


(71) 


(72) 


The  TDOAs  for  all  combinations  of  ordered  pairs  of  bearings  are  solved 
using  (71)  for  each  reflector  to  give  the  complete  set  of  possible  expected  TDOAs  which 
can  be  shown  to  number 


2  N 


M(A  +  1)] 


M 2  (A3  +2N2  +  A^)-M(A2  +  A) 


(73) 


TDOAs  where  M  is  the  number  of  transmitting  sources.  It  should  be  emphasized  that  this 
set  is  derived  using  the  DOA  estimations  from  the  first  step  of  the  proposed  scheme. 

The  set  of  expected  TDOAs  is  compared  with  the  observed  TDOAs, 
which  have  been  calculated  from  the  time  markers  using  a  weighted  version  of  the 
generalized  cross  correlation  estimator.  For  a  given  source,  the  bearing  association  is  then 
accomplished  by  matching  ordered  pairs  of  the  expected  TDOAs  to  the  observed 
TDOAs. 
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Image 


Figure  23.  Single-array,  single-source,  single-reflector  scenario  used  to  illustrate  the 
expected  TDOA  relationship  between  the  LOS  and  the  NLOS  signal. 

c.  Spatial  Bearing  Association 

As  depicted  in  Figure  24,  when  a  single  source  is  considered  and  the  exact 
footprint  of  the  reflectors  is  known,  the  association  problem  can  be  solved  spatially 
without  the  need  of  TDOA  estimation.  This  saves  processing  load  which  is  desirable  in 
wireless  sensor  network  applications.  For  the  single  reflector  scenario,  two  cases  can  be 
identified.  In  the  first,  one  bearing  is  outside  of  the  sector  defined  between  the  reflector 
and  the  array  as  in  Figure  24(a).  This  bearing  can  be  clearly  identified  as  the  LOS  signal. 
In  the  second  case,  both  bearings  reside  inside  the  sector  as  in  Figure  24(b).  In  this  case, 
the  localization  problem  (discussed  in  the  next  section)  is  solved  twice  for  each  of  the 
potential  LOS  bearings.  These  solutions  will  generate  two  distinct  points.  One  is  the 
position  of  the  source  while  the  second  is  an  image  of  the  source  related  to  the  specific 
reflector.  The  latter  point  can  be  readily  determined  to  be  unrealistic  based  on  the  known 
location  of  the  reflector  and  the  former  is  identified  as  the  actual  position  of  the  source. 
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Array  Array 


Figure  24. 


(a) 


(b) 


Spatial  bearing  association  when  one  source  transmits  and  the  exact 
footprint  of  the  reflector  is  known. 


When  multiple  reflectors  are  considered,  each  one  defines  an  independent 
sector  and,  in  the  case  where  a  sector  contains  two  bearings,  the  association  can  be  solved 
as  above.  When  each  sector  contains  at  most  one  bearing,  then  the  least  squares 
localization  problem  described  in  the  next  section  is  solved  as  many  times  as  the  number 
of  the  bearings  by  considering  a  different  bearing  as  the  LOS  each  time.  The  least  squares 
formulation,  which  converges  to  a  solution,  simultaneously  solves  both  the  association 
problem  and  the  localization  problem.  If  no  solution  results  from  the  procedure,  then  no 
LOS  signal  exists  and  the  localization  problem  is  solved  by  considering  all  bearings  as 
NLOS. 

In  the  preceding  discussion,  it  is  assumed  that  all  reflections  are  “single 
bounce.”  This  is  a  reasonable  assumption  provided  a  signal  strength  threshold  is  used  at 
each  receiver. 
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2.  Localization  using  both  LOS  and  NLOS  Signals 

The  third  step  in  the  proposed  scheme  is  the  source  localization  algorithm  itself. 
In  this  section,  we  begin  by  describing  the  single  array  case  and  then  expand  the 
localization  procedure  to  the  multiple  arrays  case.  For  the  latter,  we  propose  both 
centralized  and  distributed  solutions. 


a.  Single  Array  Localization 


As  we  have  noted  earlier,  the  proposed  scheme  is  capable  of  providing 
single  array  localization.  As  shown  in  Figure  25,  if  at  least  one  NLOS  bearing  is 
available,  the  source  position  is  found  as  the  intersection  of  the  LOS  bearing  and  the 
bearing  with  angle  26 0  -  0R  with  respect  to  the  image  array.  This  image  array  is  defined 


in  a  similar  manner  to  the  image  source  described  earlier  and  its  position  is  given  by 


7  = 

1 

_ 1 

2xra  Xa 

1 

_ 1 

2yRA-yA_ 

(74) 


We  also  note  that 


and 


yM  =  tan  0oxRA + yRo 


(75) 


tan(90-#o) 

Combining  (75)  and  (76),  we  see  that 
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tan# 
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Substituting  (77)  and  (78)  into  (74),  we  have 

(l -  tan2  0o)xA  +  2  tan 0o  (yA  - yRo ) 

_  1  +  tan2  0O 

2  tan  0oxA  +  (2  tan2  60  -l)yA-2  tan2  6oyRo 
1  +  tan2  60 


(76) 


(77) 
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(79) 


47 


Using  (42),  the  least  squares  fonnulation  of  the  single  array  localization  problem  when  a 
single  NLOS  signal  exists  in  conjunction  with  the  LOS  signal  is 


1  -  tan  0L 

1  -tan(2  0O-0R) 


LA  -  tan  6 LxK 
vA ,  -  tan  (  2  -  6 R  j  xA 


(80) 

(81) 


When  the  reflections  of  at  least  two  reflectors  are  present  the  localization 
task  can  be  performed  by  using  only  the  NLOS  signals.  Considering  the  scenario  of 
Figure  25,  matrix  A  and  the  corresponding  vector  b  in  (42)  become 


A  = 


1 

1 


-tan  (2  0ol-0R1) 
tan  ( 2  0 0 2  0 R  2 ) 


JNr  -  tan (20 o]  -0Rl)xAV 
yKV  -  tan  (2#0l  -  0RX )  xA1, 


(82) 


R1 


Figure  25.  Single-array  NLOS-only  localization  with  two  reflectors. 


b.  Proposed  Multiple-Array  Centralized  Localization  Scheme 

Consider  the  multiple  array-multiple  reflector  scenario  of  Figure  26  in 
which  K  arrays  have  available  the  NLOS  reflections  from  N  reflectors.  This  problem  is 
over-determined  and  can  be  solved  using  the  least-squares  estimator  proposed  earlier. 

The  formulation  of  the  matrix  A  and  the  corresponding  vector  b  is 
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(83) 
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where  Ak  e  C(V'lk~  is  the  matrix  of  data  obtained  by  the  kth  array  and  bk  e  C  v+I  is  the 
corresponding  vector.  These  can  be  shown  to  be 
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Figure  26.  Single-source  localization  using  multiple  arrays  in  the  presence  of 

multiple  reflectors. 

This  localization  scheme  is  a  centralized  scheme  which  involves  a  data 
fusion  center  to  process  the  incoming  data  collected  by  the  arrays.  In  a  WSN,  the  fusion 
center  may  be  an  array  of  sensors  as  well.  The  centralized  approach  has  several 
fundamental  drawbacks  when  implemented  in  WSNs.  The  most  important  is  the 
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increased  energy  consumption  of  the  fusion  center.  Thus,  the  array  assigned  to  act  as  the 
fusion  center  will  have  a  shorter  lifetime  due  to  the  added  processing  and  communication 
load  which  will  affect  the  lifetime  of  the  entire  WSN. 

c.  Proposed  Multiple-Array  Distributed  Localization  Scheme 

As  an  alternative  to  the  centralized  approach  above,  a  distributed  solution 
can  be  implemented  that  solves  the  localization  problem  for  each  single  array 
individually  and  then  averages  the  individual  array  estimates.  A  two-step  approach,  this 
solution  is  summarized  in  Figure  27.  This  distributed  approach  is  of  particular  interest  in 
wireless  sensor  network  deployments  since  it  distributes  the  processing  load  among 
multiple  sensor  nodes  and  reduces  the  communication  burden  which  is  the  primary 
energy  consumer  [2]. 


Figure  27.  Proposed  distributed  localization  scheme. 

C.  SUMMARY 

In  this  chapter,  we  proposed  a  least  squares  estimator  for  DOA -based  localization 
which  is  unbiased  when  the  noise  is  Gaussian-distributed  with  zero  mean.  This  estimator 
solves  an  over  determined  Vandennonde  system  of  equations  which  is  known  to  be 
computationally  efficient  and  accurate. 

Based  on  this  least  squares  error  estimator,  we  proposed  a  passive  source 
localization  scheme  which  exploits  the  NLOS  signals  from  non-cooperative  sources.  The 
proposed  solution  is  a  hybrid  DOA/TDOA  source  localization  scheme  and  is  comprised 
of  three  parts:  a  DOA  estimator,  an  association  algorithm  for  the  identified  signal 
bearings,  and  the  source  localization  scheme  itself.  The  recently  proposed  Space  Division 
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Multiple  Access  (SDMA)-based  receiver  was  used  for  DOA  estimation.  TDOA 
information  was  used  to  discriminate  between  the  line-of-sight  (LOS)  and  the  NLOS 
signals  and  to  associate  the  incoming  multipath  signal  with  the  corresponding  source  and 
reflector  pair.  It  was  shown  that,  in  special  cases,  the  proposed  scheme  is  capable  of 
solving  the  association  problem  spatially  without  the  need  for  TDOA  information.  A 
technique  was  also  provided  to  estimate  the  position  and  the  orientation  of  the  reflectors 
when  site-specific  database  information  is  not  available.  Both  centralized  and  distributed 
variants  of  the  proposed  scheme  were  presented  with  the  latter  being  of  particular  interest 
in  WSNs. 

In  the  next  chapter,  we  evaluate  the  performance  of  the  proposed  passive  source 
localization  scheme  using  MATLAB. 
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IV.  SIMULATION  RESULTS  OF  THE  PROPOSED 
LOCALIZATION  SCHEME 

In  this  chapter,  simulation  results  are  provided  to  validate  the  performance  of  the 
proposed  localization  scheme.  We  begin  with  a  discussion  of  the  MATLAB  simulation 
code,  the  supporting  link  budget  analysis,  and  the  performance  metrics.  Following  this, 
simulation  results  are  provided  and  analyzed  for  both  the  single  reflector  and  the  multiple 
reflector  cases. 

A.  SET-UP  OF  THE  SIMULATION 

The  performance  of  the  proposed  localization  scheme  is  validated  using  the 
MATLAB  simulation  environment.  Details  of  the  underlying  simulation  are  provided  in 
the  following  subsections. 

1.  Matlab  Simulation 

The  MATLAB  code  used  in  this  thesis  is  included  in  the  Appendix  and  is 
implemented  in  three  steps  according  to  the  block  diagram  of  Figure  19.  The  operational 
scenario  simulated  in  MATLAB  contains  one  source  and  two  reflectors  with  a 
narrowband  signal  is  transmitted  at  a  given  SNR  from  the  source.  The  propagation 
environment  is  modeled  as  discussed  in  Chapter  II  and  a  link  budget  analysis  is 
conducted.  The  incident  signal  at  the  arrays  is  processed  and  the  response  of  each  array  is 
given  by  the  function  resignal.m.  The  SDMA  receiver  is  simulated  with  the  MATLAB 
function  doasdma.m.  Once  the  DOA  estimations  have  been  detennined,  the  simulation 
associates  the  reflector-source  bearing  pairs  and  performs  the  localization.  Finally,  the 
MATLAB  code  simulates  and  compares  the  LOS-only  localization  scheme  and  the 
proposed  localization  scheme  which  utilizes  both  LOS  and  NLOS  signals. 

2.  Environment  Simulation  and  Link  Budget  Analysis 

Using  the  scenario  shown  in  Figure  28,  we  now  provide  an  example  link  budget 
analysis.  This  single-source,  single-reflector  example  can  be  expanded  to  scenarios  with 
multiple  reflectors  or  sources.  The  source  transmits  a  narrowband  signal  with  30  kHz 
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bandwidth  at  a  carrier  frequency  of  300  MHz  .  The  signal  is  received  at  the  array  via  two 
multipath  components.  The  first  path  is  the  LOS,  while  the  second  is  a  “single  bounce” 
NLOS  signal.  The  power  of  the  transmitted  signal  is  0.1  mW  .  The  source  is  located  at 
coordinates  (jcs,_ys)  =  (170  m,230  m) ,  while  the  array  is  at  the  origin.  The  reflector  has 
orientation  0X  =10°  with  respect  to  the  x  -axis.  The  reference  point  yR  is  located  at 
(0m,400m).  The  angle  0L=  45°  corresponds  to  the  DOA  of  the  LOS  path,  while 
0R  =  79.48°  is  the  DOA  of  the  reflection.  The  length  of  the  LOS  path  is  282.84  m  ,  while 
the  length  of  the  NLOS  path  is  667.08  m .  Each  individual  array  contains  50  sensors 
randomly  distributed  in  an  area  of  size  15  m  x  1 5  m . 

Assuming  a  system  temperature  of  400  K,  the  transmitted  signal  SNR  is  equal 
to  11 7.81  dB  .  The  received  SNRs  for  the  LOS  and  the  NLOS  components  are  64.93  dB 
and  49.24  dB ,  respectively.  The  NLOS  path  is  weaker  than  the  LOS  path  by  15.69  dB  . 
As  expected,  this  is  the  result  of  a  longer  propagation  path  combined  with  the  reflection 
loss. 


Figure  28.  Single-source,  single-reflector  scenario  used  in  the  reported  link  budget 

analysis. 
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3. 


Performance  Metrics 


We  model  three  uncertainties  in  the  simulation  scenarios.  The  first  is  introduced 
by  the  noise  of  the  received  signal  which  is  considered  zero-mean  Gaussian.  A  second  is 
generated  by  the  uncertainty  in  the  position  of  the  sensors  within  the  arrays.  Finally,  the 
different  arrays  are  randomly  distributed  in  the  specified  area.  To  capture  the 
performance  of  our  proposed  scheme,  we  define  a  metric  based  on  an  estimation  7  of  the 
distance  of  the  source  from  the  array  r  and  given  as 

Yr  =  'JCTI2  +  /V  (RMS  error) 

where 

|  M 

pr = — V  ( r  -  r )  (mean  error)  , 

P  i= 1 

1  M  ^ 

a2r  =  —  ^  (r  -  jur )"  (error  variance) 

P  i= 1 

and  P  is  the  number  of  Monte  Carlo  simulations. 

B.  SIMULATION  RESULTS  OF  THE  PROPOSED  LOCALIZATION 
SCHEME 

In  this  section,  simulation  results  are  provided  to  evaluate  the  perfonnance  of  the 
proposed  scheme.  In  the  first  subsection,  we  examine  the  case  of  a  single-source,  single¬ 
reflector.  In  the  second  subsection,  we  extend  this  to  multiple  reflectors. 

In  the  following  discussion,  the  term  known  reflector  implies  that  the  knowledge 
was  obtained  from  existing  databases  or  by  a  beacon,  while  the  term  unknown  reflector 
means  that  the  position  of  the  reflector  was  estimated  as  described  in  Chapter  III.  The 
mapping  of  the  reflector  is  not  free  of  errors.  A  standard  deviation  of  0.25  m  was  used  for 
the  position  of  the  known  reflector.  This  value  was  chosen  to  agree  with  the  observations 
of  [21]. 


(85) 

(86) 
(87) 


1.  Single  Source-single  Reflectors 

In  the  first  scenario,  a  single  reflector  provides  strong  multipath  components  and  a 
single  transmitting  source  with  position  (xs,ys)  =  (170  m,230  m)  exists  in  the  field.  The 
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reflector  has  orientation  6]  =10°  with  respect  to  the  x-axis.  The  reference  point  of  the 
reflector  is  yR,  located  at  (0,400  m) .  The  WSN  covers  a  square  area  of  45  mx45  m 

centered  at  the  origin.  The  general  layout  of  scenario  1  is  shown  in  Figure  29.  Seventy 
five  Monte  Carlo  simulations  were  conducted  and  a  comparison  was  made  between  the 
proposed  scheme  which  exploits  the  NLOS  signals  and  the  LOS-only-based  scheme. 
Each  array  was  comprised  of  50  sensor  elements  randomly  distributed  in  an  area  of 
225  m2 . 


Figure  29.  Scenario  1:  Single  source  located  at  (170  m,230  m) .  Single  reflector  with 
orientation  9X  =10°  and  reference  point  at  yRl  =  400  m  .  Variable  number 
of  arrays  randomly  distributed  in  a  2025  m2  area. 

The  proposed  localization  scheme  outperforms  the  LOS-only  based  localization 
as  shown  in  Figure  30.  This  is  because  more  bearings  are  available  to  the  proposed 
scheme  to  fonnulate  the  position  estimate.  Additionally,  the  condition  number  of  matrix 
A  was  found  to  be  smaller  for  the  proposed  scheme  as  shown  in  Table  4.  The  signals 
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from  a  common  source  (i.e.,  actual  source  or  reflector)  will  tend  to  be  clustered  in  a 
realistic  wireless  sensor  network  scenario  since  the  relative  distance  to  the  source  of 
interest  is  typically  large  when  compared  with  the  spacing  between  neighboring  arrays. 
However,  the  angle  separation  between  those  clusters  is  also  relatively  large.  This  is  why 
the  least  squares  problem  for  the  proposed  scheme  is  well  conditioned  compared  with  the 
LOS-only  localization  scheme 


Figure  30.  Scenario  1 :  RMS  error  for  both  the  proposed  scheme  (LOS  and  NLOS 
signals)  and  the  LOS-only  based  localization  scheme  in  the  presence  of 

known  reflectors. 


Array 

3 

4 

5 

6 

7 

HA) 

57.8158 

30.3956 

26.7982 

25.3431 

24.4903 

k(A2) 

2.5343 

2.543 

1.7376 

1.6736 

1.7588 

Table  4.  Scenario  1 :  Condition  numbers  as  a  function  of  the  number  of  arrays  for 

the  LOS-only  localization  and  the  proposed  localization  scheme  ( d2 ) , 

respectively. 
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The  effect  of  the  conditioning  of  the  least  squares  problem  for  both  the  LOS-only 
and  the  proposed  LOS-NLOS  localization  schemes  is  shown  in  Figure  3 1  where  the  RMS 
error  is  presented  as  a  function  of  the  distance  between  the  source  and  the  WSN.  Three 
arrays  of  30  sensors  each  were  used.  The  proposed  LOS-NLOS  scheme  provides  accurate 
source  position  estimates  when  the  distance  is  large  while  the  LOS-only  scheme  is 
numerically  unstable  and  inaccurate.  We  found  that  in  many  of  the  Monte  Carlo 
simulations  the  LS  problem  of  the  LOS-only  scheme  is  ill  conditioned  for  source-array 
distances  larger  than  600  meters. 


_ i _ i _ i _ i _ i _ i _ i _ i _ i _ 

300  400  500  600  700  800  900  1000  1100  1200 


Distance  (meters) 

Figure  3 1 .  Scenario  1 :  RMS  error  as  a  function  of  the  distance  of  the  source  from  the 
arrays  for  both  the  proposed  scheme  (LOS  and  NLOS  signals)  and  the 
LOS-only  based  localization  scheme  in  the  presence  of  known  reflectors 

A  comparison  between  the  centralized  and  the  distributed  localization  scheme  is 
shown  in  Figure  32.  Both  implementations  perfonn  similarly  well,  but  in  this  particular 
example,  the  distributed  approach  is  seen  to  be  slightly  more  accurate.  Additionally,  for 
the  reasons  discussed  earlier,  the  distributed  approach  is  the  more  suitable  approach  for  a 
WSN. 
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Figure  32.  Scenario  1 :  RMS  error  for  both  the  centralized  and  the  distributed 

configuration  of  the  proposed 


The  accuracy  of  the  proposed  localization  scheme  is  reduced  when  the  reflector 
position  and  orientation  is  unknown.  Figure  33.  compares  the  accuracy  of  the  proposed 
scheme  for  both  known  and  unknown  reflectors  with  the  LOS-only  based  localization 
scheme.  An  iterative  approach,  which  updates  the  estimated  reflector  position,  improves 
the  accuracy  of  the  proposed  scheme.  Figure  34  presents  this  improvement  for  the 
proposed  localization  scheme  as  a  function  of  the  number  of  iterations  when  the 
sequential  LS  approach  is  used  for  five  arrays. 
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Figure  33. 


Figure  34. 


Scenario  1 :  RMS  error  for  the  proposed  localization  scheme  using  known 
and  unknown  reflectors  and  the  LOS-only  based  localization  scheme. 


Scenario  1 :  RMS  error  as  a  function  of  the  number  of  iterations  for  the 
proposed  localization  scheme  with  an  unknown  reflector  using  the 
sequential  LS  approach  with  five  arrays. 
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As  discussed  earlier,  matrix  A  of  the  least  squares  problem  includes  errors.  Thus, 
the  proposed  localization  scheme  was  also  implemented  using  the  TLS  solution.  This 
TLS  solution  outperforms  the  LS  solution  as  shown  in  Figure  35.  This  is  in  contrast  with 
what  was  found  for  the  LOS-only  based  localization  scheme  in  Chapter  III  and  is  the 
result  of  the  well-conditioned  nature  of  the  proposed  scheme.  However,  the  improved 
location  accuracy  is  achieved  at  the  expense  of  the  larger  processing  load  required  by  the 
singular  value  decomposition  which  makes  it  difficult  to  deploy  in  a  WSN. 


Figure  35. 


Scenario  1 :  RMS  error  for  the  proposed  localization  scheme  using  least 
squares  solution  and  the  total  least  squares  solution. 


2.  Single  Source-multiple  Reflectors 

The  second  scenario  examines  the  accuracy  of  the  proposed  scheme  when  two 
reflectors  are  present.  The  geometric  configuration  is  shown  in  Figure  36  and  is  the  same 
as  in  scenario  1,  except  that  an  additional  reflector  exists  with  orientation  02  =10°  and  at 
reference  point  yR  =  -300  m  .  Again  75  Monte  Carlo  simulations  were  conducted  and  a 

comparison  was  made  between  the  proposed  scheme  which  exploits  the  NLOS  signals 
and  the  LOS-only-based  scheme.  Each  array  was  comprised  of  50  sensor  elements 
randomly  distributed  in  an  area  of  225  m2 . 
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Figure  36.  Scenario  2:  Single  source  located  at  (170  m,230  m) .  Two  reflectors,  both 

with  orientation  6  -  10°  and  reference  point 
at  yRl  =  400  m  and  yR2  =  -300  m  respectively.  A  variable  number  of  arrays 

were  randomly  distributed  in  a  2025  m2  square  area. 

The  proposed  localization  scheme  outperforms  the  LOS-only  based  localization 
as  shown  in  Figure  37  when  the  reflectors  are  known.  Furthermore,  when  analyzing  a 
LOS-based  scheme,  one  must  consider  the  case  in  which  the  scheme  incorrectly  identifies 
a  NLOS  signal  to  be  LOS.  As  shown  in  Figure  38,  this  misidentification  causes  the  LOS- 
only  approach  to  completely  break  down  and  return  erroneous  results.  In  contrast,  the 
proposed  scheme,  by  design,  effectively  associates  the  incoming  signals. 
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Figure  37.  Scenario  2:  RMS  error  for  both  the  proposed  scheme  (LOS  and  NLOS 
signals)  and  the  LOS-only  based  localization  scheme  in  the  presence  of 
one  and  two  known  reflectors 


Number  of  Arrays 

Figure  38.  Scenario  2:  RMS  error  for  both  the  proposed  scheme  (LOS  and  NLOS 
signals)  and  the  LOS-only  based  localization  scheme  when  the  latter  takes 
one  NLOS  signal  as  the  LOS  signal. 
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The  proposed  localization  scheme  can  also  be  used  for  NLOS-only  localization. 
As  shown  in  Figure  39,  these  estimates  are  significantly  improved  when  more  than  one 
reflector  is  present  (as  is  typically  the  case  in  a  realistic  scenario).  This  is  also  due  to  the 
clustered  nature  of  the  bearings,  as  discussed  earlier,  which  affect  the  conditioning  of  the 
least  squares  problem.  The  condition  numbers  of  the  A  matrix  as  a  function  of  the 
number  of  arrays  for  NLOS-only  localization  in  the  presence  of  one  and  2  reflectors  are 
shown  in  Table  5. 
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Figure  39.  Scenario  2:  RMS  error  for  the  NLOS-only  localization  of  the  proposed 
localization  scheme  in  the  presence  of  one  and  two  known  reflectors. 


Array 

3 

4 

5 

6 

7 

k(A,) 

147.7123 

98.8117 

84.1275 

76.5002 

75.3193 

*(4) 

3.7182 

3.7283 

3.7319 

3.7366 

3.7435 

Table  5.  Condition  numbers  as  a  function  of  the  number  of  arrays  for  the  NLOS- 

only  localization  of  the  proposed  localization  scheme  in  the  presence  of  one 
reflector  (A, )  and  two  reflectors ( A2 ) ,  respectively. 
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C.  SUMMARY 

In  this  chapter,  simulation  results  were  provided  to  validate  the  performance  of 
the  proposed  localization  scheme.  We  began  with  a  discussion  of  the  MATLAB 
simulation  code  and  the  performance  metrics.  Following  this,  simulation  results  were 
provided  and  analyzed  for  both  the  single  reflector  and  the  multiple  reflector  cases. 


65 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


66 


V.  CONCLUSION 


This  thesis  proposed  a  least  squares  error  estimator.  Based  on  the  proposed 
estimator,  this  thesis  developed  a  passive  source  localization  scheme  which  exploits  the 
NLOS  signals  from  non-cooperative  sources.  The  proposed  solution  is  a  hybrid 
DOA/TDOA  source  localization  scheme  and  is  comprised  of  three  parts:  a  DOA 
estimator,  an  association  algorithm  for  the  identified  signal  bearings,  and  the  source 
localization  scheme  itself.  This  scheme  allows  a  wireless  sensor  network  to  (1)  perform 
single  array  localization,  (2)  perform  the  localization  in  a  distributed  fashion,  (3)  obtain 
source  location  estimates  with  NLOS  signals  only  and  (4)  improve  the  location  estimates 
compared  with  those  obtained  using  the  LOS  information  only.  The  proposed  solution 
requires  the  knowledge  of  the  position  of  potential  reflectors  which  can  be  obtained 
dynamically  using  the  included  reflector  mapping  algorithm.  A  reflector  knowledge 
updating  procedure  based  on  the  sequential  least  squares  was  used  when  the  reflectors 
were  unknown  and  the  accuracy  of  the  proposed  scheme  was  further  improved. 
Simulation  results  were  provided  to  demonstrate  the  effectiveness  of  the  proposed 
scheme  and  the  conditioning  of  the  least  squares  problem  associated  with  it  and  compare 
it  to  existing  solutions. 

A.  SIGNIFICANT  CONTRIBUTIONS 

There  are  two  significant  contributions  in  this  thesis.  The  first  is  a  least  squares 
error  estimator  for  DOA  localization  which  is  unbiased  when  the  noise  is  Gaussian- 
distributed  with  zero  mean.  This  estimator  solves  an  over  determined  Vandermonde 
system  of  equations  which  is  known  to  be  computationally  efficient  and  accurate. 

The  second  significant  contribution  of  this  thesis  is  a  passive  source  localization 
scheme  which  exploits  the  NLOS  signals  from  non-cooperative  sources.  The  simulation 
results  showed  that  the  proposed  localization  scheme  provided  high  accuracy  estimates 
and  outperformed  the  LOS-only  based  localization  schemes.  This  is  because  more 
bearings  are  available  and  the  conditioning  of  the  least  squares  problem  is  better  for  the 
proposed  scheme 
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B.  FUTURE  WORK 

Throughout  this  thesis,  it  has  been  assumed  that  the  positions  of  the  sensors  are 
known  without  errors.  The  presence  of  errors  affects  the  accuracy  of  the  DO  A 
estimations  and  in  turn  the  accuracy  of  the  proposed  localization  scheme.  A  future  effort 
may  examine  methods  to  reduce  the  impact  of  those  errors. 

This  thesis  did  not  address  initial  signal  of  interest  (SOI)  detection.  Since  the  SOI 
SNR  at  the  receiving  sensor  network  is  dependent  to  a  large  degree  on  the  characteristics 
of  the  non-cooperative  emitter,  future  work  could  include  a  study  of  SOI  detection,  and 
subsequent  localization,  as  a  function  of  received  SNR. 

This  thesis  assumed  stationary  sources  and  arrays.  This  assumption  can  be 
extended  to  include  platforms  that  are  moving  at  low  speeds.  A  future  effort  could  also 
examine  the  perfonnance  of  the  proposed  scheme  in  tracking  applications. 

A  signal  power  threshold  was  used  to  ensure  that  all  reflections  are  “single 
bounce.”  Future  work  may  solve  the  association  problem  of  multiple-hop  reflections  and 
exploit  the  information  contained  in  those  reflections. 

Finally,  an  analysis  which  compares  the  computational  cost  of  both  the 
centralized  and  the  distributed  configurations  of  the  proposed  localization  scheme  may  be 
a  subject  for  a  future  work. 
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APPENDIX 


This  is  the  Matlab  code  used  in  this  thesis  to  evaluate  the  performance  of  the 
proposed  source  localization  scheme  in  the  presence  of  one  or  two  reflectors. 


%THESIS  MATLAB  CODE 

%Created  by  Georgios  Tsivgoulis 


clear;  clc;  elf  reset; 
%%  TRANSMITTED  SIGNAL 


f=3E8 ; 

TB=1 / (4*f ) ; 

Nchips=64 ; 

SPC  =  16; 
cp=3E8 ; 

lambda  =  cp/f; 
k  =  2*pi/lambda; 
tend=SPC*Nchips ; 
t= [0 : (tend-1) ] *TB/ (tend-1) ; 
s=sin (2*pi*f*t) ; 

S=zeros ( 1 , length (t) )  ; 

S ( 1 ,  :)  =S  ( 1 ,  : ) +exp  (lj  *  (s) ) ; 


freq  of  baseband  modulation 
time  bandwidth  product  =  4 
Number  of  Chips  /  Time  seq 
Samples  Per  Chip 
speed  of  light 
c/ f 

wavenumber 

time  snapshot  end 

discritization  of  signal 

incident  signal 


%Signal  SNR 
P=  1 0  A  -  4  ; 

K=1.38*10A-23; 

Tk=400; 

B=3E4 ; 

SNR=10*logl0 (P/ (K*Tk*B) ) 


%  signal  power  (watt) 
%  Boltzman  constant 
%  system  temperature 
%  signal  bandwidth 
%  signal  SNR 


%%  TARGET  POSITION-REFLECTOR  SIMULATION 


S-9'9'9'5'9'5'9'9'9'5'9'5'9'5'9'9'9'9'9'5'9'2'9'9'9'5'9'2'9'5'9'9'9'9'9'5'9'9'9' 

oooooooooooooooooooooooooooooooooooooooo 


%Simulation  of  the  source  and  the  arrays 

9'9'9'9'2'9'5'9'2'9'5'9'2'9'9'S-9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'2'S-9' 

ooooooooooooooooooooooooooooooooooooooooo 


tic 

Xt=170;  Yt=230; 
stdev=0 .25; 
deviation 
T= [Xt  Yt] ; 

Rt=sqrt (T (1) A2+T  (2) A2)  ; 
thol=10*pi/180; 
YRol=400; 
tho2=-5*pi/180; 
YRo2=-300 ; 


targets  position 
reflector  position 

source  distance 
reflector  1  orientation 
reflector  1  reference  point 
reflector  2  orientation 
reflector  2  reference  point 
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c=input ( ' give  the  number  of  clusters-receivers : ' ) 
array  size=input ( ' give  the  size  of  the  array:') 

N=input ( ' give  the  number  of  sensors/cluster:') 
nruns=input ( ' give  the  number  of  Monte  carlo  estimates:') 

%  array  position 
Xs=zeros ( 1 , c) ;Ys=zeros ( 1 ,  c)  ; 

%f or  i=l: (c-1) /2; 

%Ys (2*i) =15*i; 

%Xs (2*i+l)=15*i; 

%end 

for  run=l:nruns;  %  Monte  carlo  simulation  begin  from  here 

fprintf ( ' run : %2 . Of \n ' , run) ; 

YRol=YRol+stdev*randn; 

YRo2=YRo2+stdev*randn; 
for  i=l:c; 

Ys (i) =45*rand; 

Xs (i) =45*rand; 

tan_thL= (Yt-Ys (i) ) / (Xt-Xs (i) ) ; 
thL=atan(tan  thL) ; 


9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9-9'9'9'9-9'9' 

oooooooooooooooooooooooooooooo 

%Simulation  of  the  reflector 

9'9'5'9'9'9'5'9'2'9'5'9'5'9'9'9'2'9'5'9'9'9'9'9'9'9'5'9'2'S- 

oooooooooooooooooooooooooooooo 


%REFLECT0R1 

XRsl= (1/ (1+tan (thol) A2)  ) *Xs (i)  +  (tan (thol) / (1+tan (thol) A2))*(Ys(i)- 
YRol )  ; 

YRsl=tan (thol) *XRsl+YRol; 

Xsiml=2*XRsl-Xs (i)  ; 

Ysiml=2*YRsl-Ys (i)  ; 

thAl=2*thol+atan ( (Ysiml-T (2) ) / (T (1) -Xsiml) ) ;  %DOA  REF  actual  #1 
XR1 1=50 ;  YR1 l=tan (thol) *XRll+YRol; 

XR12=200;  YR12=tan (thol) *XR12+YRol; 

thl l=atan ( (YRll-Ys (i) ) / (XRll-Xs (i) ) ) +5*pi/180; 

thl2=atan ( (YR12-Ys (i) ) / (XR12-Xs (i) ) ) -3*pi/180; 

%REFLECTOR2 

XRs2= (1/ (1+tan (tho2) A2) ) *Xs (i) + (tan (tho2) / (1+tan (tho2) A2))*(Ys(i)- 
YRo2 ) ; 

YRs2=tan (tho2) *XRs2+YRo2; 

Xsim2=2*XRs2-Xs (i)  ; 

Ysim2=2*YRs2-Ys (i)  ; 

thA2=2*tho2+atan ( (Ysim2-T (2 ) ) / (T ( 1 ) -Xsim2 ) ) ;  %DOA  REF  actual  #2 
XR21=50;  YR2 l=tan (tho2) *XR21+YRo2; 

XR22=200;  YR22=tan (tho2) *XR22+YRo2; 
th21=atan ( (YR21-Ys (i) ) / (XR21-Xs (i) ) ) -5*pi/180; 
th22=atan ( (YR22-Ys (i) ) / (XR22-Xs (i) ) ) +3*pi/180; 

%%  RECEIVED  SIGNAL 


S-9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9-9'9'9'9'9' 

ooooooooooooooooooooooooooo 

%Signal  strength  simulation 

9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9-9'9'9'9'9' 

ooooooooooooooooooooooooooooo 
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%Array  Gain 
GA=10*logl0 (N) +2; 


Assume  dipole  antenna  elements 


%Free  space  path  loss 

dL=sqrt  (  (Xt-Xs  (i)  )  '“'2+  (Yt-Ys  (i)  )  A2)  ; 

LSL=-20*logl0 (lambda) +20*logl0 (dL) +21 . 98; 
dAl=sqrt ( (Xt-Xsiml) A2+ (Yt-Ysiml) A2) ; 

LSAl=-20*logl0 (lambda) +20*logl0 (dAl) +21 . 98; 
dA2=sqrt ( (Xt-Xsim2) A2+ (Yt-Ysim2) A2)  ; 

LSA2=-20*logl0 (lambda) +20*logl0 (dA2) +21 . 98; 

%Reflection  loss 

thTl=asin (sin (pi/2-thAl) / sqrt (5) ) ; 

Hrefl= (cos (pi/2-thAl) -sqrt (5) *cos (thTl) ) / (cos (pi/2- 
thAl) +sqrt (5) *cos (thTl) )  ; 

Lref l=20*logl0 (abs (Hrefl) )  ; 

thT2=asin (sin (pi/2+thA2) / sqrt (5) ) ; 

Href 2= (cos (pi/ 2+thA2 ) - 

sqrt (5) *cos (thT2) ) / (cos (pi/2+thA2) +sqrt (5) *cos (thT2) ) ; 
Lref2=20*logl0 (abs (Href2) ) ; 


%Received  SNR 
SNRL=SNR-LSL+GA; 
SNRA1=SNR-LSA1+Lref 1+GA; 
SNRA2=SNR-LSA2+Lref 2+GA; 


9'9'9'S-9'9'9'&-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9' 

ooooooooooooooooooooooooooooooooo 

%Received  Signal  simulation 

oooooooooooooooooooooooooooooooooo 


[B1,B2,B3,B]  =resignal (thL, thAl , thA2 , N, array  size,k); 
XL=B1*S;  %Received  signal 
XL=awgn (XL, SNRL, 'measured' ) ; 

XA1=B2*S;  %Received  signal 
XAl=awgn (XA1 , SNRA1 , ' measured ' ) ; 

XA2=B3*S;  %Received  signal 
XA2=awgn (XA2 , SNRA2 , ' measured ' ) ; 

X=XA1 +XA2  +XL ; 

%%  DOA  ESTIMATION 
thESTsdma=doa  sdma(X,B,N); 

%thESTsdma2=doa  sdma (XA1 , B, N) ; 

%thESTsdma3=doa  sdma (XA2 , B, N) ; 

%thESTsdma= [ thESTsdmal ; thESTsdma2 ; thESTsdma3 ] ; 

%%  BEARING  ASSOCIATION 

%LOS=f ind (thESTsdma>thll  |  thESTsdma<thl2  | thESTsdma>th2 1  | 

thESTsdma>th22  ) ; 

REFl=f ind (thESTsdma<thll  &  thESTsdma>thl2 ) ; 

REF2=find (thESTsdma>th21  &  thESTsdma<th22 ) ; 

LOS=find (thESTsdma<thl2&thESTsdma>th22) ; 
thLest (i) =thESTsdma (LOS) ; 
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thAestl ( i ) =thESTsdma (REF1 ) ; 
thAest2 ( i ) =thESTsdma (REF2 ) ; 
end 


%%  LOCALIZATION 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S' 

ooooooooooooooooooo 

%Matrix  Formulation 

S-9'9'9'5'9'5'9'9'9'2'9'2'9'9'9'9'9'2' 

ooooooooooooooooooo 

for  j=3:c; 

for  i=l  :  j 

A(i,:)  =  [l  -tan (thLest  ( i ))]  ; 

A(i+j,:)=[l  -tan (2*thol-thAestl (i) ) ] ; 
A(i+2*j,:)=[l  -tan (2*tho2-thAest2 (i) ) ] ; 
b(i, l)=Ys (i) -tan (thLest (i) ) *Xs (i) ; 
b (i+j , 1) =  Ysiml-tan (2*thol-thAestl (i) ) *Xsiml; 
b (i+2*j  ,  1) =  Ysim2-tan (2*tho2-thAest2 (i) ) *Xsim2; 

%Matrix  formulation  for  image  points 
Arefl (i, : ) = [1  -tan (thAestl (i) ) ] ; 
brefl (i, 1) =Ys (i) -tan (thAestl  (i) ) *Xs  (i)  ; 

Aref2  (i,  : )  =  [1  -tan (thAest2 (i) ) ] ; 
bref2 (i,  1) =Ys (i) -tan (thAest2  (i) ) *Xs  (i)  ; 
end 


9'9'9'9'2'9'5'9'2'9'9'9'2'9'5'9'9'9'9'S-5'9'9'9'2'9'2'9'9'S-9'S-9'9'5'9' 

oooooooooooooooooooooooooooooooooooo 

%l)LOS-only  localization 

S-9'5'9'9'9'9'9'9'9'5'9'2'9'9'9'9'9'5'9'9'9'9'9'2'9'9'9'2'9'5'9'2'9'5' 

ooooooooooooooooooooooooooooooooooo 


Xestl  =A ( 1 : j , : ) \b ( 1 : j , 1 ) 

Rtestl=sqrt (Xestl (1) A2+Xestl (2) A2) ; 

Yl=A (1 : j , : ) *Xestl ; 
kappal=cond (A (1 : j , : ) ) ; 

thetal=asin (norm (b ( 1 : j , 1 ) -Y1 ) /norm (b ( 1 : j  ,  1 ) ) ) ; 
etal=norm (A ( 1 : j , : ) ) *norm (Xestl ) / norm ( Y1 ) ; 

Condbxl (run, j -2 ) =kappal / (etal*cos ( thetal ) )  ; 

CondAxl ( run, j -2 ) =kappal+ ( kappa 1 A2*tan (thetal ) ) /etal ; 
Yt_sdmal (run, j -2 ) =Xestl ( 1 )  ; 

Xt_sdmal (run, j -2 ) =Xestl (2 ) ; 

%Bias  of  the  Target  location 
biasYt_sdmal ( run,  j -2) =T (2) -Xestl (1) ; 
biasXt_sdmal ( run, j -2 ) =T (1) -Xestl (2) ; 
biasRt  sdmal (run, j -2 ) =Rt-Rtestl ; 

% REFLECTOR  DETERMINATION% 

%REFLECT0R1 

Xestref 1  =Aref 1 ( 1 : j , : ) \bref 1 ( 1 : j , 1 ) ; 

tan_thMl= (Xestrefl (1) -Xestl (1) ) / (Xestl (2) -Xestrefl (2) ) 
tho_estl= ( (pi/2) -atan (tan_thMl) ) ; 
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tho_estl *  180/pi ; 

YRo_estl= (Xestl ( 1 ) +Xestref 1  ( 1 ) ) /2- 

tan (tho_estl) * ( (Xestl (2) +Xestrefl (2) ) / 2 ) ; 

Xestref2  =Aref2\bref2 ; 

%REFLECT0R2 

tan_thM2= (Xestref2 (1) -Xestl (1) ) / (Xestl (2) -Xestref2 (2) ) ; 

tho_est2= ( (pi/2 ) -atan (tan_thM2 ) ) ; 

tho_est2*180/pi; 

YRo_est2= (Xestl ( 1 ) +Xestref2 ( 1 )  )  / 2- 

tan (tho  est2) * ( (Xestl (2) +Xestref2 (2 ) ) / 2 ) ; 


XRsestl= (1/ (1+tan (tho_estl) A2) ) *Xs (i) + (tan (tho_estl) / (1+tan (tho_estl) A2 
)  )  * (Ys (i) -YRo_estl) ; 

YRsestl=tan (tho  estl )  *XRsestl+YRo  estl; 

Xsimestl=2*XRsestl-Xs (i) ; 

Ysimestl=2*YRsestl-Ys (i) ; 

XRsest2= (1/ (1+tan (tho_est2) A2) ) *Xs (i) + (tan (tho_est2) / (1+tan (tho_est2) A2 
)  )  * (Ys (i) -YRo_est2) ; 

YRsest2=tan (tho  est2  )  *XRsest2+YRo  est2; 

Xsimest2=2*XRsest2-Xs (i) ; 

Ysimest2=2*YRsest2-Ys (i) ; 

%Matrix  Formulation 
for  i=l : j 

Ar(i,:)=[l  -tan (thLest (i) ) ] ; 

Ar (i  +  j  ,  : )  =  [1  -tan (2*tho_estl-thAestl (i) ) ] ; 

Ar (i+2* j  ,  : )  =  [1  -tan (2*tho_est2-thAest2 (i) ) ] ; 
br (i, 1) =Ys (i) -tan (thLest (i) ) *Xs (i) ; 

br(i+j,l)=  Ysimestl-tan (2*tho  estl-thAestl (i) ) *Xsimestl; 
br(i+2*j,l)=  Ysimest2-tan (2*tho  est2-thAest2 (i) ) *Xsimest2; 
end 


S-9'9'9'2'9'5'9'2'9'9' 

ooooooooooo 

9'9'9'9'9'9'9'9-9'S-9' 

ooooooooooo 

%SCENARIO  1 


S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooo 

oooooooooooo 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%2) Scenario  1 : Centralized  proposed  localization ( 1  known  reflector) 

9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9-9'9-9'9'9'S-9'9-9'9'9'9'9'9-9'9'9'9-9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A2=A (1 : 2* j , : ) ;  b2=b ( 1 : 2* j ,  : )  ; 

Xest2=A2 \b2 
Y2=A2  *Xest2  ; 

Rtest2=sqrt (Xest2 (1) A2+Xest2 (2) A2) ; 

Yt_sdma2 (run, j -2 ) =Xest2 ( 1 )  ; 

Xt_sdma2 (run, j -2 ) =Xest2 (2 ) ; 

kappa2=cond (A2 ) ; 

theta2=asin (norm (b2-Y2 ) / norm (b2 ) ) ; 
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eta2=norm (A2 ) *norm (Xest2 ) / norm ( Y2 ) ; 

Condbx2 (run, j -2 ) =kappa2 / (eta2*cos ( theta2 )  )  ; 

CondAx2 (run, j -2 ) =kappa2+ ( kappa2  A2  *tan ( theta2 ) ) /eta2 ; 

%Bias  of  the  Target  location 
biasYt_sdma2 (run, j-2) =T (2) -Xest2 (1) ; 
biasXt_sdma2 ( run, j-2) =T (1) -Xest2 (2)  ; 
biasRt  sdma2 ( run, j -2 ) =Rt-Rtest2 ; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9-9'9'9'9'9'9'9'9-9'S-9'9'9'9'9'9'9'S-9'9'9'9'9'9-9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%3) Scenario  1 : Distributed  propossed  localization  (1  known  reflector) 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9-9'9'9'9-9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

for  i=l : j 

A3= [A (i,  : ) ; A (i  +  j ,  :  )  ]  ; 
b3= [b (i) ; b  (i  +  j )  ]  ; 

Test3 ( : , i) =A3\b3; 
end 

Xest3=mean (Test3, 2) ; 

Yt_sdma3 (run, j -2 ) =Xest3 ( 1 )  ; 

Xt_sdma3 (run, j -2 ) =Xest3 (2 ) ; 

Rtest3=sqrt (Xest3 (1) A2+Xest3 (2)  A2)  ; 

%Bias  of  the  Target  location 
biasYt_sdma3 (run, j -2) =T (2 ) -Xest3 ( 1) ; 
biasXt_sdma3 (run, j -2 ) =T ( 1 ) -Xest3 (2); 
biasRt_sdma3 ( run, j -2 ) =Rt-Rtest3 ; 

S-S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S-S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%4) Scenario  1 : Centralized  proposed  localization  (1  unknown  reflector) 

9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'S-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A4=Ar ( 1 : 2* j ,  :)  ; 
b4=br  ( 1 :  2* j ,  : )  ; 

Xest4  =A4\b4;  Y4=A4*Xest4; 

Rtest4=sqrt (Xest4 (1) A2+Xest4 (2) A2) ; 

Yt_sdma4 (run, j -2 ) =Xest4 ( 1 ) ; 

Xt_sdma4 (run, j -2 ) =Xest4 (2 ) ; 


%Bias  of  the  Target  location 
biasYt_sdma4 (run, j -2) =T (2 ) -Xest4 ( 1 ) ; 
biasXt_sdma4 (run, j -2) =T (1) -Xest4  (2) ; 
biasRt  sdma4 ( run, j -2 ) =Rt-Rtest4 ; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%5) Scenario  l:NLOS-only  localization  (1  known  reflector) 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9-9'9'9'9-9'S-9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'&-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A5=A ( j  +1 : 2* j ,  :)  ; 
b5=b ( j  +1 : 2* j , 1)  ; 

Xest5  =A5\b5; 

Rtest5=sqrt (Xest5 (1) A2+Xest5 (2) A2) ; 

Y5=A5*Xest5 ; 
kappa5=cond (A5 ) ; 

theta5=asin (norm (b5-Y5 ) /norm (b5 ) ) ; 
eta5=norm (A5 ) *norm(Xest5) /norm(Y5) ; 
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Condbx5 (run, j -2 ) =kappa5/ (eta5*cos ( theta5 ) )  ; 

CondAx5 ( run, j -2 ) =kappa5+ ( kappa5 A2*tan ( theta5 ) ) /eta5 ; 
Yt_sdma5 ( run, j -2 ) =Xest5 ( 1 )  ; 

Xt__sdma5  (run,  j  -2  )  =Xest5  (2  )  ; 

%Bias  of  the  Target  location 
biasYt_sdma5 (run, j -2) =T (2 ) -Xest5  (1) ; 
biasXt_sdma5 (run, j -2) =T (1) -Xest5 (2); 
biasRt_sdma5 ( run, j -2 ) =Rt-Rtest5 ; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'S-9'9-9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%6) Scenario  1 : Centrilized  proposed  localization  (TLS,  1  reflector) 

9'9'9'9'9'9'9'9-9'S-9'9'9'9-9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

nl  =  size (A4 , 2 ) ; 

C  =  [A4  b4 ] ; 

[U1  SI  VI]  =  svd(C, 0) ; 

Vll  =  VI ( 1 : nl , 1+nl : end) ; 

V12  =  VI ( 1+nl : end, 1+nl : end) ; 

Test6  =  -VI 1 /VI 2; 

Rtest6=sqrt (Test6 (1) A2+Test6 (2) A2) ; 

Yt_sdma6 ( run,  j -2 ) =Test6 ( 1 )  ; 

Xt_sdma6 (run, j -2 ) =Test6 (2  )  ; 

%Bias  of  the  Target  location 
biasXt_sdma6 (run, j-2)=T(2)-Test6(l)  ; 
biasYt_sdma6 (run, j-2)=T(l)-Test6(2)  ; 
biasRt_sdma6 (run, j -2 ) =Rt-Rtest6; 


S-9'9'9'9'9'5'9'9'9'9'9'2'9'9'9'9'9'9'9'9'9'9'9'5'9'5'9'2'9'9'9'9'9'5'9'5'9'5'9'5'9'9'9'2'9'9'9'9'9'5' 

ooooooooooooooooooooooooooooooooooooooooooooooooooo 

%7 ) Scenario2 : Centralized  proposed  localization  (2  known  reflector) 

9-9'9'9'9'9'9'9-9'9-9'S-9'9'9'9'9'9-9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooo 

Xest7=A\b 
Y7=A*Xest7 ; 

Rtest7=sqrt (Xest7 (1) A2+Xest7 (2) A2) ; 

Yt_sdma7 (run, j -2 ) =Xest7 ( 1 ) ; 

Xt_sdma7 ( run, j -2 ) =Xest7 (2  )  ; 

kappa7=cond (A) ; 

theta7=asin (norm (b-Y7 ) / norm (b) ) ; 
eta7=norm(A) *norm(Xest7) /norm(Y7) ; 

Condbx7 (run, j -2 ) =kappa7 / (eta7*cos (theta7 ) )  ; 

CondAx7 ( run, j -2 ) =kappa7+ ( kappa7  A2*tan ( theta7 ) ) /eta7 ; 

%Bias  of  the  Target  location 
biasYt_sdma7 ( run, j -2 ) =T (2 ) -Xest7 (1) ; 
biasXt_sdma7 (run, j -2) =T (1) -Xest7  (2) ; 
biasRt  sdma7 (run, j -2 ) =Rt-Rtest7 ; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%8) Scenario  2 : Centralized  proposed  localization  (2  unknown  reflector) 
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9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'S-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

Xest8  =Ar\br; 

Rtest8=sqrt (Xest8 (1) A2+Xest8 (2) A2) ; 

Yt_sdma8 (run,  j -2 ) =Xest8 ( 1 )  ; 

Xt_sdma8 (run, j -2 ) =Xest8 (2 )  ; 


%Bias  of  the  Target  location 
biasYt_sdma8 (run, j -2) =T (2) -Xest8 (1) ; 
biasXt_sdma8 ( run, j-2)=T(l)-Xest8(2)  ; 
biasRt_sdma8 (run, j -2 ) =Rt-Rtest8 ; 

9'9'9'9'9'9'9'9'2'9'5'9'2'9'5'9'9'9'9'9'2'9'2'9'2'9'9'9'2'9'9'9'9'9'5'9'2'9'5'9'2'9'5'9'9'9'5'9'2'9'9'9'2'9'5'9'2'9'9'9'9'9'5'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%9) Scenario  2 : Distributed  proposed  localization  (2  known  reflector) 

9'9'5'9'2'9'9'9'9'S-9'9'9'9'9'9'9'S-2'9'2'9'9'9'2'S-5'9'9'9'9'9'9'9'5'9'5'9'9'9'2'9'5'9'9'9'9'9'9'9'2'9'9'S-9'9'2'9'9'9'9'9'9'9'2' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

for  i=l:j 

A9= [A (i,  : ) ; A (i  +  j ,  : ) ; A (i+2* j ,  :  )  ]  ; 
b9= [b (i) ; b (i  +  j ) ;b (i+2* j ) ]  ; 

Test9 ( : , i) =A9\b9; 
end 

Xest9=mean (Test9, 2) ; 

Yt_sdma9 ( run, j -2 ) =Xest9 ( 1 )  ; 

Xt_sdma9 ( run, j -2 ) =Xest9 (2 )  ; 

Rtest9=sqrt (Xest9 (1) A2+Xest9 (2) A2) ; 

%Bias  of  the  Target  location 
biasYt_sdma9 ( run, j-2)=T(2)-Xest9(l)  ; 
biasXt_sdma9 ( run, j-2)=T(l)-Xest9(2)  ; 
biasRt_sdma9 (run, j -2 ) =Rt-Rtest9; 

9'9'2'9'5'9'9'9'9'9'2'9'9'9'9'9'2'9'9'9'5'9'9'9'2'9'2'9'9'9'9'9'2'9'9'9'9'9'9'S-2'9'2'S-2'9'9'S-2'9'9'9'9'9'9'9'2'9'2'9'5'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

% 1 0 ) Scenario  2:NL0S-only  (2  known  reflector) 

9'9'9'9'2'9'2'9'9'9'9'9'5'9'5'9'9'9'9'9'5'9'9'9'2'9'9'S-2'9'9'9'9'9'5'9'9'9'5'9'9'9'2'9'9'S-9'9'9'9'9'9'2'S-2'9'9'9'9'9'9'9'2'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A10=A ( j  +1 : 3* j  ,  :)  ; 
blO=b ( j+1 : 3* j , 1) ; 

XestlO  =A10\bl0; 

RtestlO=sqrt (XestlO (1) A2+Xestl0 (2) A2)  ; 

Y10=A10*Xestl 0 ; 
kappalO=cond (A10) ; 

thetalO=asin (norm (blO-YlO ) / norm (blO ) ) ; 
etalO=norm (A10) *norm (XestlO) /norm(YlO) ; 

Condbxl 0 (run, j -2 ) =kappalO/ (etalO*cos ( thetal 0 ) ) ; 

CondAxl 0 (run, j -2 ) =kappalO+ ( kappa 10 A2  *tan (thetal 0 ) ) /etalO ; 

Yt_sdmal0 (run, j -2 ) =XestlO ( 1 )  ; 

Xt_sdmal0 ( run, j -2 ) =XestlO (2  )  ; 

%Bias  of  the  Target  location 
biasYt_sdmalO (run, j-2)=T (2) -XestlO (1)  ; 
biasXt_sdmalO ( run, j-2)=T(l) -XestlO (2)  ; 
biasRt_sdmalO ( run, j -2 ) =Rt-RtestlO ; 

S-9'9'9'5'9'9'9'9'9'9'9'9'9'9'9'9'9'5'9'9'9'9'9'2'9'5'9'9'9'9'9'2'9'2'9'2'9'9'9'9'9'5'9'9'9'9'9'9'9'2'9'9'9'9'9'2'9'9'9'2'9'2'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%ll)LOS-only  localization  (fault  angle) 

9-9-9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9-9'9'9'9-9'9-9'9-9'9'9'9'9'9-9'9'9'9-9'9-9'9-9'9-9'9-9'9-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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for  i=l:j— 1; 

All  ( i+1 , : )  =  [ 1  -tan (thLest (i+1) ) ] ; 
bll(i  +  l,:)  =  [l  -tan (thLest  (i  +  1) ) ]  ; 
end 

All(l,:)=[l  -tan (thAestl (1) ) ] ; 
bll (!,:)=[!  -tan (thAestl (1) )] ; 

Xestll  =All\bll; 

Rtestll=sqrt (Xestll (1) A2+Xestll (2) A2) ; 
Yt_sdmal 1 ( run, j -2 ) =Xestl 1 ( 1 )  ; 

Xt  sdmal 1 (run, j -2 ) =Xestll (2 ) ; 

%Bias  of  the  Target  location 
biasYt_sdmal 1 ( run, j-2) =T (2) -Xestll (1) ; 
biasXt_sdmall (run, j -2 ) =T ( 1) -Xestll (2) ; 
biasRt  sdmal 1 ( run, j -2 ) =Rt-Rtestl 1  ; 

end 

end 


%%  ERROR  ANALYSIS 

%1 ) Statistics  for  LOS-only  localization 

bias  sdmal= [biasXt  sdmal'  biasYt  sdmal'  biasRt  sdmal']; 
bias 

meanbiasRt  sdmal=mean (biasRt  sdmal); 
position  error 

varbiasRt  sdmal=var (biasRt  sdmal); 
of  position  error 

rmsbiasRt  sdmal=sqrt (meanbiasRt  sdmal . A2+varbiasRt  sdmal); 
position  error 

meanCondbxl=mean (Condbxl) 
meanCondAxl=mean (CondAxl) 

%2 ) Statistics  for  Centralized  proposed  localization ( 1  known 
bias  sdma2= [biasXt  sdma2 '  biasYt  sdma2 '  biasRt  sdma2 ' ] ; 
bias 

meanbiasRt  sdma2=mean (biasRt  sdma2); 
position  error 

varbiasRt  sdma2=var (biasRt  sdma2); 
of  position  error 

rmsbiasRt  sdma2=sqrt (meanbiasRt  sdma2 . A2+varbiasRt  sdma2); 
position  error 

meanCondbx2=mean (Condbx2) 
meanCondAx2=mean (CondAx2) 


%3 ) Statistics  for  Distributed  proposed  localization ( 1  known 
bias  sdma3= [biasXt  sdma3 '  biasYt  sdma3 '  biasRt  sdma3 ' ] ; 
bias 

meanbiasRt  sdma3=mean (biasRt  sdma3) ; 
position  error 


%Matrix  of 
%Mean 
%Variance 
%RMS 

reflector) 
%Matrix  of 

%Mean 

%Variance 

%RMS 


reflector) 
%Matrix  of 

%Mean 
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^Variance 


varbiasRt  sdma3=var (biasRt  sdma3) ; 
of  position  error 
rmsbiasRt  sdma3=sqrt (meanbiasRt  sdma3 . A2+varbiasRt  sdma3);  %RMS 
position  error 


%4 ) Statistics  for  Centralized  proposed 
reflector) 

bias  sdma4= [biasXt  sdma4 '  biasYt  sdma4 ' 
bias 

meanbiasRt  sdma4=mean (biasRt  sdma4); 
position  error 

varbiasRt  sdma4=var (biasRt  sdma4); 
of  position  error 

rmsbiasRt  sdma4=sqrt (meanbiasRt  sdma4 . A 
position  error 


localization (1  unknown 
biasRt  sdma4 ' ] ;  %Matrix  of 

%Mean 

%Variance 

2+varbiasRt  sdma4);  %RMS 


%5 ) Statistics  for  NLOS-only  localization  (1  known  reflector) 
bias  sdma5= [biasXt  sdma5 '  biasYt  sdma5 '  biasRt  sdma5 ' ] ; 
bias 

meanbiasRt  sdma5=mean (biasRt  sdma5) ; 
position  error 

varbiasRt  sdma5=var (biasRt  sdma5) ; 
of  position  error 

rmsbiasRt  sdma5=sqrt (meanbiasRt  sdma5 . A2+varbiasRt  sdma5); 
position  error 


%Matrix  of 
%Mean 
%Variance 
%RMS 


meanCondbx5=mean (Condbx5) 
meanCondAx5=mean (CondAx5) 


%6) Statistics  for  Centrilized  proposed  localization  (TLS,  1 
bias  sdma6= [biasXt  sdma6'  biasYt  sdma6'  biasRt  sdma6']; 
bias 

meanbiasRt  sdma6=mean (biasRt  sdma6) ; 
position  error 

varbiasRt  sdma6=var (biasRt  sdma6) ; 
of  position  error 

rmsbiasRt  sdma6=sqrt (meanbiasRt  sdma6 . A2+varbiasRt  sdma6) ; 
position  error 


reflector) 
%Matrix  of 

%Mean 

%Variance 

%RMS 


%7 ) Statistics  for  Centralized  proposed  localization (2  known 
bias  sdma7= [biasXt  sdma7 '  biasYt  sdma7 '  biasRt  sdma7 ' ] ; 
bias 

meanbiasRt  sdma7=mean (biasRt  sdma7); 
position  error 

varbiasRt  sdma7=var (biasRt  sdma7); 
of  position  error 

rmsbiasRt  sdma7=sqrt (meanbiasRt  sdma7 . A2+varbiasRt  sdma7); 
position  error 


reflector) 
%Matrix  of 

%Mean 

%Variance 

%RMS 


meanCondbx7=mean (Condbx7) 
meanCondAx7=mean (CondAx7) 


%8 ) Statistics  for  Centralized  proposed  localization  (2  known  reflector) 
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^Matrix  of 


bias  sdma8= [biasXt  sdma8 '  biasYt  sdma8 '  biasRt  sdma8 ' ] ; 
bias 

meanbiasRt  sdma8=mean (biasRt  sdma8);  %Mean 

position  error 

varbiasRt  sdma8=var (biasRt  sdma8);  ^Variance 

of  position  error 

rmsbiasRt  sdma8=sqrt (meanbiasRt  sdma8 . A2+varbiasRt  sdma8);  %RMS 
position  error 


%9) Statistics  for  Centralized  proposed 
reflector) 

bias_sdma9= [biasXt_sdma9 '  biasYt_sdma9 ' 
bias 

meanbiasRt  sdma9=mean (biasRt  sdma9) ; 
position  error 

varbiasRt  sdma9=var (biasRt  sdma9) ; 
of  position  error 

rmsbiasRt  sdma9=sqrt (meanbiasRt  sdma9.A 
position  error 


localization (2  unknown 
biasRt  sdma9'];  %Matrix  of 

%Mean 

%Variance 

2+varbiasRt  sdma9);  %RMS 


%10) Statistics  for  NLOS-only  localization  (1  known  reflector) 

bias  sdmalO= [biasXt  sdmalO'  biasYt  sdmalO'  biasRt  sdmalO'];  %Matrix  of 

bias 

meanbiasRt  sdmalO=mean (biasRt  sdmalO);  %Mean 

position  error 

varbiasRt  sdmalO=var (biasRt  sdmalO);  ^Variance 

of  position  error 

rmsbiasRt  sdmalO=sqrt (meanbiasRt  sdmalO . A2+varbiasRt  sdmalO) ;%RMS 
position  error 


meanCondbxl 0=mean (CondbxlO) 
meanCondAxl 0=mean (CondAxlO) 


%1 1 ) Statistics  for  LOS-only  localization  (fault  angle) 

bias  sdmall= [biasXt  sdmall'  biasYt  sdmall'  biasRt  sdmall'];  %Matrix  of 
bias 

meanbiasRt  sdmall=mean (biasRt  sdmall);  %Mean 

position  error 

varbiasRt  sdmal l=var (biasRt  sdmall);  %Variance 

of  position  error 

rmsbiasRt  sdmal l=sqrt (meanbiasRt  sdmal 1 . A2+varbiasRt  sdmall) ;%RMS 
position  error 


toe 


SN=3 : 7 ; 
figure ( 1 ) 

plot ( SN, rmsbiasRt  sdmal b--A SN, rmsbiasRt  sdma2,'r-*') 
xlabel ( ' Number  Arrays ' , ' Font size ' , 12 ) 
ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend (' LOS-only  localization Proposed  localization  (known 
reflector)  '  ) 
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figure  (2 ) 

plot ( SN, rmsbiasRt  sdma2  ,' g--A ', SN, rmsbiasRt  sdma3,'r-*') 
xlabel (' Number  of  Arrays Fontsize 12 ) 
ylabel ( ' \gamma-RMS (meters) ' , ' Fontsize ' , 12) 

legend (' Centralized  proposed  localization Distributed  proposed 
localization  '  ) 

figure  ( 3 ) 

plot ( SN, rmsbiasRt  sdmal b-A SN, rmsbiasRt  sdma2,'g- 

o ', SN, rmsbiasRt  sdma4, ' r-*') 

xlabel (' Number  of  Arrays Fontsize ', 12 ) 

ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend (' LOS-only  localization Proposed  localization  (known 
reflector) ' , . . . 

'Proposed  localization  (unknown  reflector) ') 
figure ( 4 ) 

plot ( SN, rmsbiasRt  sdma4 ,' r--A ', SN, rmsbiasRt  sdma6,'r-*') 
xlabel (' Number  of  Arrays Fontsize ', 12 ) 
ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend ( ' LS  proposed  scheme (unknown  reflector)  ' ,  '  TLS  proposed 
scheme (unknown  reflector) ' ) 

figure ( 5 ) 

plot ( SN, rmsbiasRt  sdmal ,' b--A ', SN, rmsbiasRt  sdma2,'r- 
SN, rmsbiasRt  sdma7, ' r-*') 
xlabel ( ' Number  Arrays ' , ' Fontsize ' , 12 ) 
ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend (' LOS-only  localization Proposed  localization  (1  known 
reflector) ' , . . . 

'Proposed  localization  (2  known  reflector) ') 
figure (6) 

plot ( SN,  rmsbiasRt  sdmall,  'b--A ',  SN,  rmsbiasRt  sdma8,'r- 

* ' , SN, rmsbiasRt  sdma4 , ' g--o ' ) 

xlabel (' Number  of  Arrays ',' Fontsize ', 12 ) 

ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend (' LOS-only  localization  (fault  angle) ',' Proposed  localization  (1 
unknown  reflector)  '  ,  .  . . 

'Proposed  localization  (2  unknown  reflector) ') 

figure ( 7 ) 

plot ( SN, rmsbiasRt  sdma5 ,' b--A ', SN, rmsbiasRt  sdmalO, ' r-* ' ) 
xlabel ( ' Number  Arrays ' , ' Fontsize ' , 12 ) 
ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend (' NLOS-only  localization  (1  reflector) ', 'NLOS-only  localization 
(2  reflector) ' ) 
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This  is  the  Matlab  code  used  in  this  thesis  to  evaluate  the  performance  of  the 
proposed  source  localization  scheme  in  the  presence  one  unknown  reflector  using  the 
sequential  least  squares. 


%MODEL  WHICH  APPLIES  THE  SEQUENTIAL  LEAST  SQUARE  IN  THE  LOCALIZATION 


clear;  clc;  elf  reset; 


Var=0.25*piA2/180A2; 
f=3E8 ; 
modulation 
TB=1 / (4*f ) ; 

=  4 

Nchips=64 ; 
seq 

SPC  =  16; 
cp=3E8 ; 

lambda  =  cp/f; 
k  =  2*pi/lambda; 

tend=SPC*Nchips ; 

t= [0 : (tend-1) ] *TB/ (tend-1) ; 

signal 

s=sin (2*pi*f*t) ; 

S= zeros ( 1 , length (t) ) ; 

S ( 1 ,  :)  =S  ( 1 ,  : ) +exp  ( 1 j  *  ( s ) )  ; 

%Signal  Power 

P=10A-4; 

K=1.38*10A-23; 

Tk=400; 

B=3E4 ; 

SNR=10*logl0 (P/ (K*Tk*B) ) 


%  freq  of  baseband 
%  time  bandwidth  product 
%  Number  of  Chips  /  Time 
%  Samples  Per  Chip 

%  c/f 

%  wavenumber 

%  time  snapshot  end 
%  discritization  of 

%incident  signal 


%%  TARGET  POSITION-REFLECTOR  SIMULATION 
tic 

Xt=170;  Yt=230;  %  targets 

position 

stdev=0 .25; 

T= [Xt  Yt] ; 

Rt=sqrt (T (1) A2+T  (2) A2)  ; 
thol=10*pi/180; 

YRol=400; 

c=input ( ' give  the  number  of  clusters-receivers : ' ) 
array  size=input ( ' give  the  size  of  the  array:') 

N=input ( ' give  the  number  of  sensors/cluster:') 
nruns=input ( ' give  the  number  of  Monte  carlo  estimates:') 
it=input ( ' give  the  number  of  iterations') 

Xs=zeros (l,c) ;Ys=zeros (l,c)  ; 
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Monte  carlo  simulation  begin  from 


for  run=l:nruns; 
here 

fprintf ( ' run : %2 . Of \n ' , run) ; 

YRol=YRol+stdev*randn; 

%Array  positions 
for  i=l:c; 

Ys (i) =45*rand; 

Xs (i) =45*rand; 

tan_thL= (Yt-Ys (i) ) / (Xt-Xs (i) ) ; 
thL=atan(tan  thL) ; 

S-S-9-9'9-9'9-9'9'S-9'9'9'9'9-9'9'9-9-9'9'9'9'9-9-9'9' 
ooooooooooooooooooooooooooo 

%Simulation  of  the  reflector 

9'9'2'9'9'9'9'9'2'9'5'9'9'9'9'9'2'9'5'9'9'9'2'9'9'9'9'9'9'S- 
oooooooooooooooooooooooooooooo 

XRsl= (1/ (1+tan (thol) A2)  )  *Xs (i)  +  (tan (thol) / (1+tan (thol) A2))*(Ys(i)~ 
YRol )  ; 

YRsl=tan (thol) *XRsl+YRol; 

Xsiml=2*XRsl-Xs (i)  ; 

Ysiml=2*YRsl-Ys (i)  ; 

thAl=2*thol+atan ( (Ysiml-T (2) ) / (T (1) -Xsiml) ) ;  %DOA  REF  actual  #1 
XR1 1=50 ;  YR1 l=tan (thol) *XRll+YRol; 

XR12=200;  YR12=tan (thol) *XRl2+YRol; 

thl l=atan ( (YRll-Ys (i) ) / (XRll-Xs (i) ) ) +5*pi/180; 

thl2=atan ( (YR12-Ys (i) ) / (XR12-Xs (i) ) ) -3*pi/180; 

S-S-9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooo 

%Signal  strength  simulation 

9'9'9'9'9'9'9'S-9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9-9'9'9' 

ooooooooooooooooooooooooooooo 

%Array  Gain 

GA=10*logl0 (N) +2 ;  %Assume  dipole  antenna  elements 

%Free  space  path  loss 

dL=sqrt ( (Xt-Xs (i) ) A2+ (Yt-Ys (i) ) A2) ; 

LSL=-20*logl0 (lambda) +20*logl0 (dL) +21 . 98; 
dAl=sqrt ( (Xt-Xsiml) A2+ (Yt-Ysiml) A2) ; 

LSAl=-20*logl0 (lambda) +20*logl0 (dAl) +21 . 98; 

%Reflection  loss 

thTl=asin (sin (pi/2-thAl) / sqrt (5) ) ; 

Hrefl= (cos (pi/2-thAl) -sqrt (5) *cos (thTl) ) / (cos (pi/2- 
thAl) +sqrt (5) *cos (thTl) ) ; 

Lref l=20*logl0 (abs (Hrefl) ) ; 

%Received  SNR 
SNRL=SNR-LSL+GA; 

SNRA1=SNR-LSA1+Lref 1+GA; 


9-9-9'9'9'9'9'9'9'9-9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9' 

ooooooooooooooooooooooooooooooooo 

%Received  Signal-  DOA  Estimation 

9'9'9'9'9'9'9'9'9'9'5'9'9'9'9'9'2'9'5'S-9'S-5'9'2'9'9'9'2'S-9'9'2'9' 

oooooooooooooooooooooooooooooooooo 

for  j=l:it; 

[B1,B2,B]  =resignal (thL, thAl , N, array  size,k); 
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XL=B1*S;  %Received  signal 
XL=awgn (XL, SNRL, 'measured'); 
XA1=B2*S;  %Received  signal 
XAl=awgn (XA1 , SNRA1 ,  ' measured ' ) ; 
X=XA1+XL; 

thESTsdma=doa  sdma(X,B,N); 


S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

oooooooooooooooooooooooo 

%Assosiation 

9-9'9'9'9'9-9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooo 

REFl=find (thESTsdma<thll  &  thESTsdma>thl2 ) ; 

LOS=f ind (thESTsdma<thl2 )  ; 
thLest (i, j ) =thESTsdma (LOS) ; 
thAestl ( i , j ) =thESTsdma (REF1 )  ; 
end 
end 

S-S-S-S-S-S-S'S-S'S-S-S-S'S-S-S-S'S-S- 

ooooooooooooooooooo 

%Matrix  Formulation  for  LOS  and  REF 

S-9-9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooo 

for  j=l:it; 
for  i=l:c; 

Alos(i+c*j-c,:)=[l  -tan (thLest (i, j )  )  ]  ; 

bios (i+c* j -c,  1) =Ys (i) -tan (thLest (i , j ) ) *Xs  (i ) ; 

Aref (i+c*j-c,:)=[l  -tan (thAestl (i, j ) ) ] ; 

bref (i+c*j-c, 1) =Ys (i) -tan (thAestl ( i , j ) ) *Xs (i) ; 

end 

end 

Xestlos ( : , 1 ) =Alos (l:c,:)\blos(l:c, :); 

Xestref ( : , 1 ) =Aref ( 1 : c,  : ) \bref ( 1 : c,  : )  ; 

9'9'2'9'2'9'9'9'2'9'9'9'2'9'5'9'2'9'2'9'2'9'2'9'2'9'2'9'2'9'2'9'2'9'9'S-9'9'5'S-2'9'9'9'2'9'2'S-2'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Reflector  estimation  using  sequential  least  square 

ooooooooooooooooooooooooooooooooooooooooooooooooooo 

for  j  =1 : c* ( it-1 ) 

C  (1 :  j+c-1, 1 : j+c-1)  =  (1/Var) *eye(j+c-l)  ; 

Sigmalos (:, : ) =inv (Alos (1 : j +c-l, :) '*C(l:j+c-l,l:j+c-l) *Alos (1 : j+c- 
1,  :)  )  ; 

Klos ( : , j ) = (Sigmalos ( : , : ) *Alos ( j+c, : ) ' ) / (Var+Alos ( j+c, : ) * 

Sigmalos ( : ,  : ) *Alos (j+c, : )  ' )  ; 

Xestlos  ( : , j+1) =Xestlos ( : , j ) +Klos (:,j) * (bios (j+c, 1)- 
Alos  (j+c, : ) *Xestlos ( : ,  j )  )  ; 

Sigmaref ( : , : ) =inv (Aref (1 : j+c-1, :) '*C(l:j+c-l,l:j+c-l) *Aref (1 : j+c-1, :) ) ; 

Kref ( : , j ) = (Sigmaref ( : , : ) *Aref (j+c, : ) ' ) / (Var+Aref (j+c, : ) * 

Sigmaref ( : ,  : ) *Aref (j+c, : )  '  )  ; 

Xestref ( : , j+1) =Xestref ( : , j ) +Kref ( : ,  j ) * (bref (j+c,  1)  - 
Aref (j+c, : ) *Xestref ( : ,  j  )  )  ; 

tan_thM= (Xestref (1,  j ) -Xestlos (1,  j) ) / (Xestlos (2,  j) -Xestref (2,  j) )  ; 
tho_est ( j ) = ( (pi/2) -atan (tan_thM) ) ; 

YRo_est ( j )  =  (Xestlos (1, j ) +Xestref (1,  j )) 72- 

tan  (tho_est ( j ) ) * ( (Xestlos (2, j ) +Xestref (2,  j ) ) / 2 )  ; 

for  i=l:c; 

XRsest ( j , i) = (1/ (1+tan (tho_est ( j ) ) A2) ) *Xs (i) + (tan (tho_est ( j ) ) / (1+tan (tho 
_est  ( j ) ) A2) ) * (Ys (i) -YRo_est ( j ) ) ; 
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YRsest ( j , i) =tan (tho_est ( j ) ) *XRsest ( j ) +YRo_est (j  )  ; 

Xsimest ( j , i) =2*XRsest ( j  ,  i) -Xs (i)  ; 

Ysimest ( j , i) =2*YRsest ( j , i) -Ys (i)  ; 

A r ( i  ,  : )  - [  1  -tan (thLest  (i) ) ]  ; 

Ar (i+c, : ) = [1  -tan (2*tho_est ( j ) -thAestl (i) ) ] ; 
br (i, 1) =Ys (i) -tan (thLest (i) ) *Xs (i) ; 

br (i+c, 1) =  Ysimest ( j , i) -tan (2*tho_est ( j ) -thAestl (i) ) *Xsimest ( j , i) ; 
end 

Xest2=Ar\br 

Rtest2=sqrt (Xest2 (1) A2+Xest2 (2) A2) ; 
biasRt  sdma2 (run, j ) =Rt-Rtest2; 
end 

%Los  localization 
Xestl=Xestlos ( : ,  1) 

Rtestl=sqrt (Xestl (1) A2+Xestl (2) A2) ; 
biasRt  sdmal (run) =Rt-Rtestl ; 
end 

meanbiasRt  sdmal=mean (biasRt  sdmal); 
position  error 

varbiasRt  sdmal=var (biasRt  sdmal); 
of  position  error 

rmsbiasRt  sdmal=sqrt (meanbiasRt  sdmal . A2+varbiasRt  sdmal); 
position  error 

meanbiasRt  sdma2=mean (biasRt  sdma2); 
position  error 

varbiasRt  sdma2=var (biasRt  sdma2); 
of  position  error 

rmsbiasRt  sdma2=sqrt (meanbiasRt  sdma2 . A2+varbiasRt  sdma2); 
position  error 

rmsbiasRt  sdmal=rmsbiasRt  sdmal *ones ( 1 , 8 ) 

SN=1 : 8; 
figure  ( 1 ) 

plot ( SN, rmsbiasRt  sdmal ,' b--A ', SN, rmsbiasRt  sdma2,'r-*') 
xlabel (' Number  of  iterations Fontsize ', 12 ) 
ylabel ( ' \gamma, RMS (meters) ' , ' Fontsize ' , 12) 

legend (' LOS-only  localization LS  proposed  (known  reflector) ') 


%Mean 

%Variance 

%RMS 


%Mean 

%Variance 

%RMS 
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This  is  the  function  used  to  determine  the  received  signal  in  the  array 

function  [B1,B2,B]  =resignal (thLOS , thREFl , N, array  size,k  ) 


Ntheta=91;  Nphi=3601;  %  Set  up  angle  matrices 

theta  =  linspace ( 0 , pi/2 , Ntheta) ; 
phi  =  linspace ( -pi , pi , Nphi ) ; 

DTheta=  theta (2) -theta (1) ;  DPhi  =  phi (2) -phi (1) ; 

[Phi,  Theta]  =  meshgrid (phi , theta) ; 

CT=cos (Theta) ;  ST=sin (Theta) ;  CP=cos(Phi);  SP=sin(Phi); 

Rx=ST.*CP;  Ry=ST.*SP;  Rz=CT;  %  unit  radial  vector 

components 


%  in  wavelengths  (lambda)  square  grid 
min  spacing  =array  size/50  ;  %  l/10th  wavelength 

sepn  min  =  0;  %  initialize  prior  to  1st  test 

xp  =  array  size  *  rand(l,N); 
yp  =  array_size  *  rand(l,N); 

% - Stats  of  Random  Array - 


while  (sepn  min<  min  spacing) 
wavelength 


min  spacing  =  . 1 


rp  =  [xp; yp] ' ; 
sepn  =  pdist(rp); 
elements  in  Random  Array 
sepn  min  =  min (sepn); 
zsepn  =  squaref orm ( sepn) ; 
array  to  square 

[Rmin  Cmin]  =  f ind ( zsepn==sepn  min,l, 
if  (sepn  min<  min  spacing) 

xp (Rmin) =array  size  *  rand; 


yp (Rmin) =array  size  *  rand; 


%  dist  between  all 

%  find  minimum  separation 
%  convert  linear  sepn 

'first');  %returns  index 
%  then  Move  one  of  them 


end 


m=f ind (theta==pi/2 )  ; 
Rxm  =  Rx  (m,  :  )  ; 

Rym  =  Ry (m, : ) ; 


want  row  for  xy  plane 
lxNphi  vector 


%  Calculate  the  actual  received  signal  for  each  angle  using  the  N 
element 

%  receiver  and  sum  all  N  signals, 
for  n=l : N 

Bll(l,n)  =  exp  (j * (k* (cos (thLOS) *xp (n) +sin (thLOS) *yp (n) ))) ; 
B22(l,n)  =exp(j*(k*(  cos (thREFl) *xp (n) +sin (thREFl) *yp (n) ))) ; 
end 


for  n=l : 3601 ; 

an(n)=(n/10-0.1) -180; 
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th (n) =an (n) *pi/180; 
for  i=l : N 

B ( i ,  n)  =exp ( 1 j  * ( k* (xp ( i ) *Rxm (n) +  yp ( i ) *Rym (n) ) ) ) ; 

end 

end 

B1=B1 1 ' ; 

B2=B22 ' ; 


This  is  the  function  used  to  estimate  the  DOA  of  the  incident  signal  using  the  SDMA 
receiver. 

function  thESTsdma=doa  sdma(X,B,N) 

%%  SDMA  RECEIVER 

Nchips=64 ; 

SPC  =  16; 

H  =  hadamard (Nchips ) ; 
codes= [ ] ; 
for  i=l : N 

codes ( i ,  : ) =H  ( i  +  1 ,  : )  ; 

end 

%  spread  codes  by  SPC 
codes  =  repmat (  codes, SPC, 1) ; 
codes  =  reshape (codes, N, SPC*Nchips) ; 
codes  =  codes*pi/2; 
codes=codes ' ; 
r=codes*X; 

v=codes*B; 

R=con j (v) '  *r ; 

R=sum (R, 2 ) ; 
maxR=max (abs  (R)  )  ; 

Psdma=abs (R) /maxR; 

Psdma= [ zeros ( 900 , 1 ) ;  Psdma ( 902 : 2702 ,:) ;  zeros ( 900 , 1 )] ; 

%DOA  SDMA 
Nphi=3601 ; 

phi  =  linspace ( -pi , pi , Nphi ) ; 
for  i=l : Nphi 
if  Psdma ( i) >0 . 5 ; 

Psdmal ( i ) =Psdma ( i ) ; 

else 

Psdmal ( i ) =0 ; 

end 
end 

POSsdma  =  imregionalmax ( Psdmal ) ; 
local  maxima->BW  array  "l"s 
PEAKsdma  =  f ind ( POSsdma) ; 

Maxima 


%  MATLAB  rtn  that  finds 
%  index  pointers  to  local 


%  Samples  Per  Chip 
%  use  Walsh  Sequence 

%  don't  want  1st  row 


%  repeat  by  spread  factor 
%  move  like  chips  together 
%  pi/2  since  inside  exponential 
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[bs , is ] =sort ( Psdma ( find ( POSsdma) ),' descend ') ;  %  sort  by  highest  peak- 
>lowest 

ThESTsdma=  PEAKsdma ( is ) ;  %  index  to  peaks 

ThESTsdma=  180*phi (ThESTsdma) /pi;  %  Angle  where  peaks  occur 

thESTsdma=ThESTsdma (1,1:2); 
thESTsdma=thESTsdma*pi/ 180; 


This  is  the  function  used  to  estimate  the  DOA  of  the  incident  signal  using  the  MUSIC 
algorithm. 

do  a  music,  m 


function  thESTmusic=doa  music (X, A, N) 


Nchips=64 ; 

SPC  =  16; 
tend=SPC*Nchips ; 

Nphi=3601 ; 

phi  =  linspace ( -pi , pi , Nphi ) ; 
Rxx=X*X ' / tend; 

[V,  Dia] =eig (Rxx) ; 
matrix 

[Y, Index] =sort (diag (Dia) ) ; 

EN=V ( : , Index ( 1 : N-2 ) ) ; 
for  i=l : 3601; 

MU (i) =abs (A ( : , i) ' *EN*EN ' *A( :  ,  i)  )  ; 
end 

Pmu=l . /MU; 
maxPmu=max (Pmu) ; 

Pmu= [ Pmu (1802:3601)  Pmu (1:1801) ] ; 
Pmu  s i c = Pmu / max  Pmu ; 

Pmusic= [zeros (1, 900)  Pmusic (901 : 27 
POSmusic  =  imregionalmax (Pmusic)  ; 
local  maxima->BW  array  "l"s 
PEAKmusic  =  find (POSmusic) ; 
local  Maxima 

[bs,is]=sort (Pmusic (find (POSmusic) 
>lowest 

ThESTmusic=  PEAKmusic ( is ) ; 
ThESTmusic=  180*phi (ThESTmusic) /pi 

occur 

thESTmusic=ThESTmusic (1,1:3); 
thESTmusic=thESTmusic*pi/180 ; 


%  Number  of  Chips  /  Time  seq 
%  Samples  Per  Chip 


%autocorrelation  matrix  of  signal 
%eigenvalues  of  autocorrelation 

%eigenvalues  in  descending  order 
%noise  subspace  matrix 


%music  pseudospectrum 


1 )  zeros (1,900) ] ; 

%MATLAB  rtn  that  finds 

%  index  pointers  to 

, 'descend');  %  sort  by  highest  peak- 

%  index  to  peaks 

;  %  Angle  where  peaks 
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